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ABSTRACT 

The upper atmospheres of close-in gas giant exoplanets ("hot Jupitcrs") are subjected to intense 
heating and tidal forces from their parent stars. The atomic (H) and ionized (H + ) hydrogen layers 
are sufficiently rarefied that magnetic pressure may dominate gas pressure for expected planetary 
magnetic field strength. We examine the structure of the magnctosphcrc using a three-dimensional 
(3D) isothermal magnetohydrodynamic model that includes: a static "dead zone" near the magnetic 
equator containing gas confined by the magnetic field; a "wind zone" outside the magnetic equator 
in which thermal pressure gradients and the magneto-centrifugal-tidal effect give rise to a transonic 
outflow; and a region near the poles where sufficiently strong tidal forces may suppress transonic 
outflow. Using dipole field geometry, we estimate the size of the dead zone to be several to tens of 
planetary radii for a range of parameters. Tides decrease the size of the dead zone, while allowing the 
gas density to increase outward where the effective gravity is outward. In the wind zone, the rapid 
decrease of density beyond the sonic point leads to smaller densities relative to the neighboring dead 
zone, which is in hydrostatic equilibrium. To understand the appropriate base conditions for the 3D 
isothermal model, we compute a simple one-dimensional (ID) thermal model in which photoelectric 
heating from the stellar Lyman continuum is balanced by collisionally-excited Lyman a cooling. This 
ID model exhibits a H layer with temperature T ~ 5, 000 — 10, 000K down to a pressure P ~ 10 — 100 
nbar. Using the 3D isothermal model, we compute maps of the H column density as well as the 
Lyman a transmission spectra for parameters appropriate to HD 209458b. Line-integrated transit 
depths ~ 5 — 10% can be achieved for the above base conditions, in agreement with the results of 
Koskinen et al. A deep, warm H layer results in a higher mass-loss rate relative to that for a more 
shallow layer, roughly in proportion to the base pressure. Strong magnetic fields have the effect of 
increasing the transit signal while decreasing the mass loss, due to higher covering fraction and density 
of the dead zone. Absorption due to bulk fluid velocity is negligible at linewidths > 100 km s^" 1 from 
line center. In our model, most of the transit signal arises from magnetically confined gas, some of 
which may be outside the LI cquipotcntial. Hence the presence of gas outside the LI equipotential does 
not directly imply mass loss. We verify a posteriori that particle mean free paths and ion-neutral drift 
are small in the region of interest in the atmosphere, and that flux freezing is a good approximation. 
We suggest that resonant scattering of Lyman a by the magnctosphere may be observable due to 
the Doppler shift from the planet's orbital motion, and may provide a complementary probe of the 
magnctosphere. Lastly, we discuss the domain of applicability for the magnetic wind model described 
in this paper as well as the Roche-lobe overflow model. 



1. INTRODUCTION 

Hot Jupiters, the gas giant exoplanets found at small 
orbital radii, D < 0.1 AU, present an opportunity to 
study planets orbiting in an extreme environment very 
close to their parent stars. They experience insolation 
levels ~ 10 4 times greater than solar system gas giants. 
As a consequence of the high temperatures generated by 
EUV heating, the gas pressure scale height is compara- 
ble to the planetary radius, creating an extended upper 
atmosphere of gas potentially observable through trans- 
mission or reflection spectroscopy. Heating may also cre- 
ate an outward pressure force, complemented by outward 
centrifugal and tidal forces, which may drive an outflow 
leading to mass and angular momentum loss from the 
planet. 

While most observational effort for the transit- 
ing planets has centered on the lower atmosphere, 
near the photospheres for optical and infrared contin- 
uum radiation at mbar-bar pressures, there has also 
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been progress in probing higher levels in the atmo- 
sphere. As these observations probe a larger area 
surrounding the planet, the transit depth is corre- 
spondingly larger than that from the photospheric ra- 
dius, R p h- Bound-bound atomic lines in the trans - 
miss ion spectrum from Nal dCharbonneau et al.l 120021) . 
HI A Vidal-Madiar et all 120031 120041: lEhrenreich et al 



12008ft. OI and CII (IVid al-Madiar et all 12004(1 . and SiHI 
(jLinskv et aIll20Toh have been observed for HD 209458b ; 
HI in HP 189733b dLecavelier Pes Etangs et al.l l2010Tl; 
and Mgll (jFossati et alj2010ft . and possibly other species 
at lower signal to noise, in WASP-12b. Bound-free tran- 
sitions from a population of hydrogen in the n = 2 
state have be en detected in the tra nsmission spectrum of 
HP 209458b (jBallester et alJl2007ft . althoug h no bound- 
bound transitions from n = 2 were found (jW inn et aLl 

12551 . 

The large transit depths observed for the Lyman a 
line in HP 209458b imply neutral H atoms occupy an 
occulting area equivalent to an optically thick disk ex- 
tending out to several planetary radii, comparable to 
the Roche-lobe radius for this close-in planet. This 
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led iVidal-Madiar et all ()2003[) to suggest that the H 
atoms are escaping from the planet, as they would no 
longer be bound to the planet outside the Roche-lobe 
radius. The presence of heavy atoms (C, O, etc.), 
whose scale height should be far smaller than that of 
the H atoms, was suggested to be due to a combina- 
tion of efficient turbulent mixing and drag forces from 
the esca ping gas, entraining the heavier elements in a n 
outflow (IVidal-Madiar et al.ll2004: iGarcia Mu5odl2007t ). 
IBen-.Taffell (j2lCT used archival HST data to rederive the 
transit depth, and challenged the claim that atmospheric 
escape was occurring, due to a derived transit depth 
that would imply the H atoms are within the Roche 
lobe radius. Two subseq uent studies dBen- Jaffell 120081 : 
IVidal-Madiar et all 120081 ) elucidated the dependence of 
transit depth on the analysis met hod and wavelength 
range studied. The recent study by iLinskv et al.l ()2010l ) 
attempts to constrain the wind velocity in HD 209458b 
using the line profiles for the CII and Silll lines. 

Most theoretical efforts to explain the aforementioned 
observations invoke EUV heating to temperatures ~ 
10 4 K, creating a thermally driven, p ossibly transonic 
hydrodynamic outflow from the planet ([Y clle 2004. 120061: 
Tian et al.ll2005HGarcia Muhozll2007tlMurrav-Clav et al.l 
2009tlStone fc Progall2009t) . While initiated by outward 
gas pressure gradients, a thermally-driven wind can be 
significantly accelerated by the tidal gravitjQ. An alter- 
native model i s to treat atmospheric escape as Roche- 
lobe overflow (|Gu et al.l 12001 iLi et al.l l2010t lLai et al.l 
assuming that the fluid can only reach the sound 
speed in the vicinity of the LI Lagrange point. 

For gas giants, H is the most abundant element and 
heating is dominat ed by phot oionization of atomic H in 
the thermosphere (Yclle 2004|). Photoionization leads to 
significantly ionized gas at high altitudes that is sub- 
ject to magnetic forces. In light of the high ionization 
fraction, the planetary magnetic field, due to dynamo 
action in the planet's core, may be dynamically impor- 
tant and affect density and velocity profiles in the upper 
atmosphere. This paper is a first attempt to include the 
effects of the planet's magnetic field combined with tidal 
and centrifugal forces on upper atmosphere structure for 
hot Jupiters. 

We model the magnetized, hot, rotating portion of the 
planet's upper atmosphere in the context of magneto- 
hydrodynamic (MHD) outflows, the theory of which was 
originally developed for stars and accretion disks. Figure 
[1] shows a cartoon of the expected structure for the wind 
from an isolated object with a surface dipolc field: the 
polar region supports an outflow while the equatorial re- 
gion contains static, magnetically confined gas (see also 
lYelld I2004D . This configuration applies interior to the 
magnetosphere-stellar wind interaction, for planets not 
too near the parent star. Near the parent star, tidal 
forces are an important consideration, and we will show 
that tides may strongly affect the density profile and size 
of the static region, and may even shut the wind off in 
the polar region. 

In the first half of the paper (§[2] through §[7]), we de- 
velop a general theory of isothermal hot Jupiter magnc- 

2 In this paper we will use the phrase "tidal gravity" to include 
the effect of both the stellar tide and the centrifugal force, assuming 
synchronized rotation (see §[5]) 




Fig. 1. — Schematic model for magnetic field structure (thick 
lines and arrows) and fluid velocity (thin arrows). In the polar 
"wind zone", magnetic field lines are open allowing outflow. In 
the equatorial "dead zone" close to the planet, the gas has zero 
velocity and no outflow occurs. The dead zone ends in a cusp- type 
neutral point, denoted by the dashed circle at several planetary 
radii, outside of which field lines are open at all angles. This figure 
is characteristic of the weak tide limit, while in the strong tide limit 
the polar wind would be partially suppressed (see §0. 

tosphcres. Section [2] discusses the problem setup and ap- 
proximations used to obtain a solution. Planetary mag- 
netic field strengths predicted by dynamo models are dis- 
cussed in section [3J Section @] discusses qualitative fea- 
tures of atmospheric structure, and motivation for the 
existence of a dead zone. Section [5] reviews the centrifu- 
gal and tidal forces, and discusses the projection of these 
forces along magnetic field lines. The structure of the 
dead zone is discussed in section [51 followed by the wind 
zone in section [7J 

The second half of the paper (§ [5] through § ll"2")) uses 
the Lyman a transmission spectra to constrain parame- 
ters of the global 3D models. In section[51 we construct a 
ID model in hydrostatic, thermal, and ionization balance 
to compute appropriate values for the base pressure and 
temperature for the 3D global models that we present in 
section [5] Mass loss rates are computed and spin-down 
torques are estimated, in section [TO] Neutral H column 
density maps for the global models are used to illustrate 
the dependence on key parameters in section 1111 Lyman 
a transmission spectra in comparison with observations 
and scattering of stellar Lyman a by the planet are pre- 
sented in section Q21 Finally, we compare and contrast 
our magnetic wind model with the standard Roche-lobe 
overflow model in section 1131 Summary and discussion 
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are given in section [TU The MHD wind equations and 
ion-neutral coupling are discussed in Appendices [3] and 
iBl respectively. 

2. PROBLEM SETUP AND APPROXIMATIONS 

In this section we outline the problem to be solved, and 
the simplifying assumptions used to find solutions. Ex- 
cept for section [51 this paper discusses a 3D isothermal 
model of the upper atmosphere. The isothermal model 
is parametrized by an effective sound speed, as well as 
the pressure at a fiducial radius. Appropriate values for 
these parameters are discussed in a simple ID spherical 
model in section [8j including photoionization and ther- 
mal equilibrium. 

In the 3D model, we compute approximate solutions 
of the one-fluid MHD equations with an inner boundary 
at the base of the warm H layer, and an outer boundary 
which extends to at least ten planetary radii (or one stel- 
lar radius) . We treat the gas as having constant isother- 
mal sound speed a = / pm p , where T is gas temper- 
ature of the fluid, [i is the mean molecular weight, fc& is 
Boltzmann's constant, and m p is the proton mass. At the 
inner boundary, the density and pressure are assumed to 
follow cquipotcntials. We assume the planet's rotation 
rate is synchronized with its orbital motion around the 
parent star, giving orbital and spin angular velocity f2. 
We work in a coordinate system centered on the planet 
and rotating at rate H. The stellar gravity is included 
in the tidal approximation, and an effective potential U 
is composed of the planetary gravity, stellar gravity, and 
the centrifugal force. We specify a specific magnetic field 
geometry which is dipole near the planet and matches 
onto a radial field at the dead zone radius. We com- 
pute the ionization fraction with a simple, optically thin 
model applied to the derived gas densities of the MHD 
model. To summarize, we have made several simplifying 
assumptions in order to focus on the new physics aris- 
ing from MHD effects. We now discuss these simplifying 
approximations in more detail. 

The simultaneous inclusion of photoionization heat- 
ing, chemical reactions and collisional coupling between 
different species, stellar tidal forces, and the simultane- 
ous interaction with the stellar wind in the presence of 
magnetic field is a formidable problem. Our approach 
is to first ignore the interaction with the stellar wind, 
but to include the effect of the stellar tidal forces felt by 
the planet's atmosphere. The interaction with the stel- 
lar wind may alter the results of this paper in several 
ways (see e.g. , iMurrav-Clav et all 120091 : iStone k Progal 
12009ft . The stellar wind will limit the size of the mag- 
netosphcre, as d etermined by stress balance at the mag- 
netopause (e.g.. IPreusse et al~ll2007| ). Reconnection be- 
tween field lines in the stellar wind and magnctosphere 
may lead to magnctosphcric convection, lim iting the hig h 
density region to be inside a plasmapause (|ParksH20"04l ). 
as for Earth. Finally, reconnection may also generate 
non-thermal plasma populations. We do not consider the 
interaction with the stellar wind in order to construct the 
simplest possible model. 

Another key approximation is that we treat photoion- 
ization heating as being spherically symmetric, creating 
a hot layer uniformly over the planet. In reality, the 
night side temperature and ionization state may depend 
on day-night heat redistribution, and downward heat 



conduction along field lines not in the planet's shadow. 
In perfect MHD, such redistribution would be highly 
constrained in the magnetically dominated upper atmo- 
sphere, but finite co nductivity may allow field lines to 
slip through the gas (|Goldll959T ). Even on the day side, 
large gas density outside the Roche lobe may project a 
non-spherically symmetric shadow on the deeper layers. 

Near the planet, the dynamo- generated, roughly dipole 
field from the planets core is expected to dominate. Mov- 
ing outward, currents generated in the magnetosphere 
comb the field lines into a nearly radial direction beyond 
the dead zone radius. Such a geometry h as been used 
before in the context of the stellar w ind (|Mestell 119681 : 
iMestel k SpruitHl987t lOkamotolH97l and we will adopt 
it here. This field geometry will be implemented in the 
global models presented in section 9, and is motivated in 
the discussion of Appendix A. 

Finally, the one-fluid approximation assumes that 
mean free paths are sufficiently small that relative motion 
of different species can be ignored. In Appcndix[B]we will 
check this assumption a posteriori for our models, which 
estimate particle densities and velocities in the dead and 
wind zones. Specifically, we will show that the electron- 
proton-hydrogen atom gas is well coupled collisionally 
for the parameters of interest, and therefore the drift ve- 
locity is small and hydrogen atoms have short mean free 
paths in the dead zone region. As a consequence, neutral 
hydrogen atoms do not fly ballistically through the mag- 
netosphere, and photoionization equilibrium is a good 
approximation when computing the ionization fraction. 
Further, we compute the rate at which magnetic field can 
drift relative to the fluid. For this thermal population of 
particles in the magnetosphere, we find that the dead 
zone gives the largest observable transit signal, and that 
the bulk of the hydrogen atoms in the dead zone are not 
escaping. 

In the next section we will review expected magnetic 
field strengths for hot Jupiters. 

3. EXPECTED MAGNETIC FIELD STRENGTHS 

The importance of the magnetic field for the upper at- 
mosphere depends critically on the field strength. How- 
ever, the magnetic fields of hot Jupiters are currently 
unconstrained by observation. This section will use the- 
oretical considerations to estimate likely field strengths 
f or hot Jupiters. 

iSanchez-Lavc ga (2004) computed that Raylcigh num- 
bers in hot Jupiters are typically much larger than 
the critical Raylcigh number for thermal convection in 
the metallic core. Using estimates of the fluid veloc- 
ity carrying the heat flux, he found that the magnetic 
Reynolds number is much larger than unity, and that 
dynamo action can occur. He argued that if the dy- 
namo operated with Elsasser number of order unity, then 
B ~ (2 J oriAs) 1 / 2 , where p is the mass density and Xb is 
the magnetic diffusivity. The dominant scaling impor- 
tant for hot Jupiters is then with rotation: B cx J7 1 / 2 . 
For synchronized planets with orbital periods of a few 
days, this scaling predicts that the field for hot Jupiters 
should be smaller than that of Jupiter (equatorial field 
-Bj,eq = 4.3 G) by a factor of a few. 

The op posite conclusion may be drawn from the recent 
results of lChristensen et al.l (|2009f ). As the rotation rate 
is increased above a critical value, the field strength no 
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longer increases with rotation rate, and dynamo simu- 
lations give a magnetic field strength B ~ (/O-^corc) 1 ^ 3 ' 
nearly independent of rotation rate and magnetic dif- 
fusivity, where F rn r P is the heat flux escaping from the 



conducting core. IChristensen et afl ([20091 ) show that 



this scaling applies to both planets and rapidly rotat- 
ing, low mass stars over many orders of magnitude in 
heat flux. They argue that the dependence on F coro 
arises since it is the heat flux reservoir that sustains the 
magnetic field against Ohmic dissipation. To be in the 
saturated regime, the Rossby number Ro must satisfy 
Ro = V e d/m £ 0.1, where V c d is the typical velocity of 
the eddies transporting heat, and £ is th e size of the con- 
ductin g region; £ ~ R for M p > 0.5Mj. iSanchez-Lavegal 
(|2004[) estimated synchronized planets in few day orbits 
to have Ro < 0.1, using heat fluxes comparable to that 
of Jupiter. Hence these planets are expected to be in the 
saturated regime. 

The large radii of hot Jupiters are currently not well 
understood, since the cooling and contraction time for 
passively cooling planets, even allowing for irradiation 
by the star, is far short er than the age for a numbe r 
of observed planets (e.g.. iFortnev fc Nettelmannl 120091 ). 
This has led to the suggestion that these planets arc 
not passively cooling, but rather have an anomalous 
source of internal heating, which is as yet unidentified 
but balances the co re cooling rate. To ass ess the re- 
quired heating rates, lArras k, Socrates! (|2009( ) computed 
cooling flux from the core for planets as a function 
of radius (their Figure 11; see lArras fc Bildstenl HM 
for a discussion of the cooling luminosity of irradiated 
hot Jupiters). For Jupiter- mass planets in the radius 
range i? p h = 1.3 — l.bRj, the cooling flux is larger 
tha n that of Jupiter by a factor 10 2 — 10 3 , for which 
the IChristensen et al.l (|2009f) scaling would give magnetic 
fields 5—10 times larger than Jupiter. Larger mass plan- 
ets with the same radius would have larger cooling fluxes, 
and vice versa. For instance, WASP 12b, WASP 17b 
and TRES 4 have radii i? p h ~ 1.8Rj and masses in the 
range 0.5 — 1.5Mj0, for which the cooling fluxes would 
be 10 3 — 10 4 times higher than Jupiter, implying fields 
larger than Jupiter by factors of 10 — 20. 

In summary, the recen t results on dynamo theory from 
IChristensen et ahl ([20091 ) . and the assumption that hot 
Jupi ter cores are subject to an externally powered heat- 
ing ([Arras fc Socrate s 2009), argue that field strengths 
may be up to an order of magnitude larger than that 
of Jupiter. We will take this as motivation to explore a 
wide range of parameter space for the magnetic field in 
our calculations. 

In the next section, we show that the photoionized 
H and H + layers are magnetically dominated for field 
strengths comparable to Jupiter or Saturn, and motivate 
the existence of a dead zone by a toy problem. 

4. DEAD ZONE- WIND ZONE STRUCTURE OF THE UPPER 
ATMOSPHERE 



lYelii ([2001 and iGarcfa Mufiol ([20071) presented de- 
tailed calculations of the transition between the molec- 
ular lower atmosphere (H 2 ), the layer dominated by 
atomic hydrogen (H), and the ionized upper atmosphere 
(H + ). In this paper, we restrict attention to the H and 

3 http:/ /exoplanet.eu/catalog-transit.php 



H + regions, where the transmission spectrum is formed. 
The strong heating in these layers due to UV photon en- 
ergy deposition raises the temperature to T ~ 10 4 K. 
As a consequence of the increased temperature and low 
mean molecular weight, the radial extent of the H and H + 
layers (> Rj) is expected to be much larger than that of 
the H2 layer above the photosphere (< (0.1 — 0.2) x Rj). 

We define the lower boundary of our wind model to be 
at the base of the warm (T > 5000 K) H layer, at base 
radius R and base pressure Pbasc- The base of the warm 
layer is a crucial parameter for the transi t depth. As dis- 
cussed in the phenomenological model of fKoskinen et ahl 
(|2010D . the transit depth of HD 209458b could be under- 
stood as being due to thermal ~ 10 4 K H gas extending 
down to ~ 10—100 nbar pressures. In section we 
compute a simple ID model for ionization and thermal 
equilibria in the H and H + layers which shows that such 
base conditions are indeed possible. 

The Lyman a transmission spectrum of HD 209458b 
shows absorption by the planetary atmosphere at 
lincwidths Av > 100 km s _1 from line center. At 
this linewidth, the cross section is er„ ~ 2 x 10~ 19 cm 2 
(see Figure [TB"]) . Optical depth unity requires a hydro- 
gen column Njj — 1/ct„ — 5 x 10 18 cm~ 2 . Assum- 
ing the gas is dominated by atomic hydrogen, the pres- 
sure at this level in the atmosphere is P ~ gm p NH — 
2 nbar (.g/300 cm s~ 2 ) (see Figure [TU|) . This is a fac- 
tor ~ 10 8 more rarefied than the optical photosphere 
for continuum radiation at pressure P p h ~ 100 mbar. 
Magnetic forces dominate in this layer if B 2 /8tt > 
2 nbar (g/300 cm s~ 2 ), implying a critical field strength 
Bcrit > 0.25 G (s/300 cm 2 s" 1 ) 1 / 2 . This is less than 
Jupiter's equatorial magnetic field Bj eq = 4.3 G and 



comparable to Saturn's equatorial field _Bs,i 



0.22 G. 



Moving upward, if the gas pressure drops much faster 
than magnetic pressure, the atmosphere can become 
highly magnetically dominated — a magnctosphere. 

The theory of thermally and magneto-centrifugally 
driven MHD winds gives guidance on the upper atmo- 
sphere structure i n the magnet ically dominated case (for 
a good review see lSpruit[|1996[ ). Consider a thought ex- 
periment in which a non-magnetic spherically symmetric 
wind with velocity Woo and mass loss rate M exists at 
time t < 0, and at time t = a dipole magnetic field 
is turned on. On which field lines can the wind over- 
power the magnetic forces and open the field lines to in- 
finity? The magnetic pressure on the equator (8 = tt/2) 
is weaker than at the footpoint (at angle 9b) by a factor 
[B{Tr/2)/B(6 b )} 2 = [R/r{ir/2)f = sin 12 6 b (see eq[T2]). 
This powerful dependence on footpoint position means 
there are always field lines near the magnetic poles which 
open to infinity on which a wind can outflow (see Fig- 
ure Q] for a cartoon) . The reason is that at the equa- 
tor, the wind ram pressure can overcome the steeply 
falling magnetic pressure at sufficiently large distance 
from the planet. If we assume the wind ram pressure 
decreases outward as pv 2 = Mvoo/A-kt 2 , then the criti- 
cal footpoint angle inside of which a "polar wind" occurs 

is sin6>b ~ {Mv oa /R 2 B 2 \ . For fiducial polar mag- 
netic field Bq = 8.6 G, R — lARj, (constant) flow speed 
«oo = 10 km s~ 4 and mass loss rate M = 10 11 g s~ 4 (mo- 
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tivated by the studies of |Yelld|200j; iGarcia Munozll2007l; 
iMurrav-Clav et aLll2009ft . we find the polar cap size oc- 
cupied by open field lines is Ob — 14°, and the last closed 
field lines is at equatorial radius r(%/2) ~ 1QR. This 
static, closed field li ne region is referred to as the "dead 
zone" ()Mestellll968l ). These estimates suggest that the 
region within a few planetary radii, where the transit sig- 
nal arises, is filled primarily with static gas in the dead 
zone, rather than outflowing g has been previously 
assumed. 

For close-in, tidally-locked planets, tides and centrifu- 
gal forces play a key role in upper atmosphere structure. 
In the next section we review the effective potential near 
the planet, and the projection of forces along dipole field 
lines. 

5. TIDAL FORCE AND MAGNETIC GEOMETRY 

Gas in the upper atmosphere of the planet is subject 
to three potential forces: gravity from the planet, tidal 
gravity from the star, and the centrifugal force due to 
the planetary rotation. For a synchronized planet, the 
centrifugal and tidal forces, or just tidal force for short, 
are comparable in strength, although their angular de- 
pendence is different. For position vector x = (r, 0, <p) 
relative to the center of the planet, and star at position 
x+ = (D,7r/2, 0) = De x , the sum of the three accelera- 
tions is 



a(x) = -VU(x), 
where the effective potential is given by 



U(x) = - 



GM P 
GM V 



GM+x ■ x+ 1 
~^7F 2 



(1) 



fl 2 r 2 (/sin 2 (9- l) 



(3) 

r 2 x ' 

Equipotentials are shown in Figure 6.2 of iKopall (|1978t ). 
The longitude-dependent function / = 1 + 3 cos 2 0. The 
vector angular velocity of the orbit is ft = £le z , where 
e z is normal to the orbital plane, and f2 = [G(M+ + 
Mpj/D 3 ] 1 / 2 . The first and second terms in eqf2]are the 
potential of the planet and star, respectively. The third 
term in cqf2] is due to the motion of the origin of the 
coordinate system. The last term in eqf2] is due to the 
centrifugal force. It may be shown that eqf2]is equivalent 
to the usual Roche potential with origin at the center of 
mass by combining the third and fourth terms. The form 
in eqf5] is an expansion in the limit r <C D, and agrees 
with t hat for "Hill's limit" found in lMurrav fc Dermottl 
(|2000t) when evaluated in the orbital plane (0 = tt/2). 
The accelerations are given by 



dU _ GM i 
dr 

— — = ffl 2 r sin 9 cos 9 
r o9 

1 dU 



+ il 2 r (/sin 2 6» - f) 



3f2 2 r sin sin cos < 



r sm ( 



(4) 
(5) 
(6) 



The vector acceleration a should not be confused with 
the sound speed a. If we denote the coordinate along 
the star-planet line x = r sin 9 cos 0, and normal to the 
orbit plane z = r cos 9, then the tidal force is outward 
when 3x 2 > z 2 , i.e. within latitudes — ir/3 to ir/3 of the 



equator. The tidal force is zero in the y-direction. Hence 
when observing a planet in the plane of the sky during 
transit, the tidal forces are inward along e z , zero along 
e y , and away from the planet along e x , the line of sight 
during transit. 

The LI and L2 Lagrangian points at radii rxi and rL2 
are found by using eqfj] along the star-planet line (6 = 
71-/2,0 = 0), giving 

rti - r L2 s r L = D(M p /U'Q 1/3 = (GMp/Sfl 2 ) 1 / 3 (7) 

The near equality t*li ~ 7"L2 is due to Mp/M* <C 1 in 
cq0 The photospheres of the observed planets arc inside 
the L1-L2 radii. 

How does the magnetic field alter the radius beyond 
which the net gravity points outward? What is needed is 
the projection a\\ = a b along field lines, where b = B/B 
is the unit vector along the magnetic field direction. In 
dipole geometry, approximately correct near the planet, 
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sm ( 



e r cos ( 



and the unit vector is 

. 1 / . sin(9 

TV [ e r COsd + e S—^~ 

where the normalization factor is 



(8) 



(9) 



N = J cos 2 9 + sin 2 0/4 = J 1 - 3 sin 2 0/4. (10) 



The parallel acceleration is then 
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cos 
N 



f GMp 
\ r 2 



n 2 r 



The quantity in brackets in ea llll must be positive in 
order for the net acceleration to be away from the planet. 
To solve for the radius at which the net acceleration is 
zero, we express in terms of r along dipole field lines 
using 



r(0)=R 



sin 2 
sin 2 0(, 



r oq sin 



(12) 



where 



r cq = i?/sin 2 6 (13) 

is the equatorial radius of the field line with footpoint at 
0£,. The "magnetic Roche lobe radius", trb, at which the 
projected acceleration an = is given by the solution of 
the equation 



GM n 



' KB 



« 2 rRB = |/« 2 ^. 



(14) 



Solutions can only exist in the range of radii 2r eq /3/ < 
fRB < ?*eq- Solutions for trb exist first at the looptop 
fRB = ^oq for loops of critical size 



' cq,crit — 



(J/2- 1/3)1/3' 



(15) 



That is, when the loop size becomes larger than about the 
L1-L2 distance, the outer part of the loop can have net 
acceleration pointing away from the planet. Note that 
this statement applies even in the plane where cos 2 = 0, 



6 



where the radial component of the tidal force points in- 
ward. The magnetic geometry allows a "magnetic Roche 
lobe" trb ~ 7*r to exist at all longitudes, in contrast to 
the unmagnetized case. 

In the next section we model the gas density in the 
dead zone. 

6. THE DEAD ZONE 

The 3D MHD wind equations are presented in Ap- 
pendix [X] There we derive the Bernoulli constant along 
field lines and discuss how gas pressure discontinuities at 
the dead zone - wind zone boundaries give rise to current 
sheets which alter the magnetic field configuration. For 
the present section which discusses the dead zone, the 
main concept needed is that hydrostatic balance applies 
along field lines. We will approximate the field lines as 
dipolar. 

In the dead zone, the velocity along field lines v = 0. 
Setting v = in eq lA2l and dotting this equation with 
b to eliminate the Lorentz force we find the equation of 
hydrostatic balance along field lines 



1 dP 
p ds 



,d\i\p 
ds 



dU 
ds 



(16) 



where d/ds = b ■ V is the derivative along field lines. 
Under the isothermal assumption, cq |16l can then be in- 
tegrated to give the run of pressure and density along a 
field line with base position (fit,, 4>) : 



P(r,6,(t>) _ P {r,e,4>) 



P{R,0 b 



p(R,0 b 



(17) 



= CXp 



fU(r,6,4>) - U(R, 6 , < 



We treat the density and pressure at the base as being 
along equipotentials. Defining p ss = P ss /a 2 = p(r = 
R 7 8 = 7r/2, <f) = 0) to be the value at the substellar point 
at the base radius, the density at the base radius at other 
points is 



pb{0b,4>)=Pss exp 



U{R, 9b,4>)-U(R,ir/2,0) 



Combining eqfTT] and [TS] then gives 
P(r,M) p(r,M) 



Pss 



(18) 



(19) 



: exp 



U{r,0,$)-U{R,ir/2,Q) 



Note that because we have assumed the density and pres- 
sure surfaces at the base are along equipotentials, e 0.1191 
satisfies 



0=-VP - pVU 



(20) 



in all three directions, not just along field lines. In this 
case —~VU is balanced by gas pressure forces due to a 
non-spherical distribution of mass. If, on the other hand, 
temperature, density or pressure surfaces were not along 
equipotentials, then eg 1201 would not be satisfied perpen- 
dicular to field lines, and the trans-field force balance 
(eq |A12[) would be required to understand the required 
currents. 
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Fig. 2. — Values of A calculated using cq |22l for the current 
database of transiting exoplanets. The parameters used are a = 
9.3 km s —1 , appropriate to T = 10 4 K, fi = 1, and Ft = -R p h- The 
open circles, X's and open triangles are planets with M p /Mj > 2, 
1/2 < M p /Mj < 2 and M p /Mj < 1/2, respectively. Values of 
^orbi JWp a nd -Rph are t aken from the current database of transiting 
exoplanets ( |http: / /exoplanet.eu/[ |. 

The potential difference in eg J19I can be written in di- 
mcnsionless form 



U(r,9,4>) - f/(.R, tt/2,0) 



(21) 



1-f 



r 



where we have defined the ratio of escape to thermal 
speed 



A: 



GM P 
Ra? 



;9.3 



0.7Mj 



l.ARj 
R 



10 km 8~ 



(22) 



and the ratio of rotation speed (or tidal potential) to 
thermal speed 



rip 

a 

: 0.043 



3.5 days 

Povb 



R 



lARj 



10 km s" 



(.23) 



The photoionization model in section [8] shows outward 
increase in T and decrease in p, implying larger a with 
radius. The density will then decrease more slowly than 
in the isothermal model for the same quantities at the 
base. 

Figures [5] and El show the values of A and e versus planet 
orbital period for the transiting exoplanets, taking M p , 
transit radius i? p h and orbital period P or b = 2ir/Q from 
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Fig. 3. — Values of e calculated using eg 1231 for the current 
database of transiting cxoplanets (http://cxoplanct.eu/). The 
open circles, X's and open triangles arc planets witn M p /Mj > 2, 
1/2 < Mp/Mj < 2 and M p /Mj < 1/2, respectively. The param- 
eters used are a = 9.3 km s _1 , appropriate to T = 10 4 K, fj, = 1, 
and R = R p ^ . 

the Extrasolar Planets Encyclopedic^. The temperature 
and mean molecular weight have been set to fiducial val- 
ues of fj, = 1 and temperature T = 10 4 K, giving sound 
speed of a = 9.3 km s _1 . 

Note that a substantial number of planets have A = 
2 — 10, implying the scale height of the gas is large enough 
that the density decrease in the dead zone is only by a 
factor of 10— 10 4 , far less than for planets with cold upper 
atmospheres more distant from their parent star. This 
increased density leads to the possibility that hot Jupiter 
upper atmospheres may be collisional to large distances 
from the planet, i.e. that the exobase, if it exists at all, 
is at radii r » R ph (|Tian et al.ll2005l iMurrav-Clav et al.l 
[2001 . 

Next, note that for the same fiducial molecular weight 
and temperature, the strength of the tide, e, is in the 
range 0.1 — 1 for a substantial number of planets. For 
large e, the typical rotational speed of a synchronized 
planet is comparable to the sound speed, or equivalcntly, 
the free fall speed in the tidal potential is comparable to 
the sound speed. Since e is evaluated at the base, the 
tidal force will dominate even more at larger distances 
from the planet. 

Figure 0] shows the run of gas and magnetic pressure 
along the equator (9 = 7r/2) as a function of radius along 
the star-planet line (</> = 0). EqfT9l was used for the 
gas pressure, and eqE] for the magnetic pressure. The 
magnetic pressure is parametrized by the plasma (3 at 
the substellar point at the base: 

4 http://exoplanct.cu/ 




Fig. 4. — Gas and magnetic pressure, normalized to the base 
gas pressure P SB = P(R,n/2,0), as a function of equatorial ra- 
dius r eq for the isothermal model. The tidal potential is evalu- 
ated along the star-planet line cos 2 (f> = 1 at the equator 9 = n/2. 
The two groups of lines starting from r = R and P = P ss are 
P(r eq , tt/2, 0)/P sb evaluated for A = 5,10. For each group, the 
line style gives the value of e. The three lines sloping down to the 
right are B 2 (r cq , n/2, 0)/(8nP as ) = (i?/r eq ) 6 //3 for the three dif- 
ferent values of equatorial /3 = 10 — 4 , 10 -2 , 1. The cusp radius in 
Figure [5] is given by the intersection of the gas and magnetic pres- 
sure curves. The gas pressure decreases outward faster for larger 
A. Beyond the Roche radius, gravity effectively points outward 
and the gas pressure begins to increase outward. For larger e, the 
Roche radius moves inward. 



The different lines in Figure [4] show gas pressure and 
magnetic pressure for different A, e, and /3 as a function 
of radius. 

For small e, the density decreases outward, and even- 
tually becomes a constant. When the tidal force is in- 
cluded the density increases outward for radii outside the 
magnetic Roche radius fea flT)) . since the sign of gravity 
points outward there. This is a dramatic effect for close- 
in planets, whose Roche radii are at only a few plan- 
etary radii, and may lead to hydrogen densities orders 
of magnitude larger than the e = case. The tidal 
gravity plays a role similar to the centrifugal force in 
models of the closed fi eld line regions in the solar wind 
(jMestel fc Spruitjll987l ) and in Jupiter's magnetosphere 
outside the corotation radius. To specify the current and 
field distributions required for this support involves a so- 
lution of the trans-field equation, which is beyond the 
scope of this work. However, such support should be 
possible when mag netic pre s sure d omin ates gas pressure. 

Nex t, we follow iMestell (| 19681 ) and IMestel fc Spruit] 
(119871) to estimate t he siz e of the dead and wind zones. 
IPneuman fc Koppl (| 1 9 7 If) showed that the dead zone 
ends at the equator in a cusp, i.e. the field in the dead 
zone approaches zero toward the cusp. Also, for v > a, 
the wind zone pressure can be neglected compared to the 
dead zone pressure, leading to the condition (also see our 
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Fig. 5. — Cusp radius as a function of equatorial plasma /3 at 
the substellar base of the atmosphere for different values of A = 
GM p /(Ra 2 ) and e = (QR/a) 2 . Here A is roughly the ratio of escape 
speed to isothermal sound speed, and e is the roughly the ratio of 
rotation speed to isothermal sound speed. Larger A (e) implies the 
density decreases outward faster (slower). The cusp radius moves 
outward (inward) for larger A (e). 



discussion leading up to eq |A14l in Appendix A) 



P 



dead 



d2 
wind 



(25) 



To determine the cusp radius, which we call for dead 
zone, we use eq[19]for the left hand side and eq]8jfor the 
right hand side. The resulting equation is 



P exp 



-All-" 

r 



3 

r 



R 2 



1 



Eg 1261 shows that the cusp radius will depend on <p m 
general due to the tidal potential. When tides can be 
ignored (e -C 1), a reasonable approximation is r^/i? ~ 
(eV/3) 1 / 6 . 

Figure [5] shows cusp radii as a function of /3 for different 
tidal strength (e) and binding parameter (A). Consider 
the fiducial case with e = 10 -3 , A = 10 and (3 = 10~ 2 . 
Initially the magnetic pressure is much larger than the 
gas pressure. The more rapid decrease of magnetic pres- 
sure implies equality at the cusp radius, r^jR — 10. In- 
creasing e, the gas pressure is larger and the cusp radius 
moves inward. The cusp radius moves outward with in- 
creasing magnetic field, i.e. decreasing (3. 

Under what conditions does a magnetosphere not 
form? Inspection of the (3 = 1 and A = 5 curves in 
Figure U shows that the magnetic and gas pressures 
arc initially equal at the base, but the gas density de- 
creases more slowly and so the magnetic field never dom- 
inates. Ignoring tides, there is an analytic criterion for 
the critical (3 CI a = /3 C rit(A) for the formation of a dead 
zone. This criterion is found by requiring simultane- 
ously P = B 2 /8tt and dP/dr = d/dr{B 2 /8tt), and yields 



P < /3crit = (6/A) 6 exp(A - 6). For A = 5, 10, 15, the val- 
ues are (3 cr a = 1.1, 2.5 and 33, agreeing with the cutoffs 
in Figure [5] Hence for large A, the gas at the base need 
not be magnetically dominated in order for the gas well 
above to base to become so. 

In summary, we have found cusp, or dead zone, radii in 
the range of a few to tens of radii R for the expected range 
of A, f3 and e. Since transit observations to date probe 
the high density gas within a few planetary radii, these 
observations may be probing static gas trapped within 
the magnetosphere, as opposed to outflowing gas. 

Nevertheless, as we will argue in section[7l a wind zone 
should exist, and we investigate its structure in the next 
section. 

7. THE WIND ZONE 

In this section wc show how the (slow magneto-) sonic 
point is affected by the magnetic geometry. Readers un- 
familiar with the MHD wind equations can consult Ap- 
pendix [A] for a brief summary of the equations. We sim- 
plify the problem by assuming that the sonic point is 
close enough to the planet for magnetic stresses to domi- 
nate over hydrodynamic stresses — the rigid field line ap- 
proximation. In this situation the fluid nearly corotates 
with the planet, and can be accelerated like "beads on a 
wire" by the magnetic field. T his "magneto -centrifugal" 
effect from stellar wind theory (jMestel 1963) becomes im- 
portant when the rotation velocity approaches th e sound 
speed at the sonic point (|Mestel fc Spruit! Il987j ). For a 
synchronized hot Jupiter, tidal forces are of comparable 
size as centrifugal forces, and the condition that the ro- 
tation velocity is comparable to the sound velocity at the 
sonic point is equivalent to the Roche-lobe radius being 
near the sonic point. The rigid field line assumption will 
typically break down near the Alfven point, which we 
estimate lies well outside the dead zone radius for most 
latitudes. 

Plugging parallel vel ocity , v = vb , and the no 
monopoles condition, eq lA71 into cq lAl| the continuity 
equation assumes the simple form 



7 (81) = 



(27) 



so that pv/B is constant on field lines, and has the inter- 
pretation of the mass loss rate per unit of magnetic flux 
alon g a fl ux tube of area oc B . Similarly, the projection 
of eq lA2l along the field can be rewritten as 



dv 
ds 



,dlnp 

ds 



dU 
ds ' 



(28) 



w hich has the Bernoulli integral along field lines (see 
oci lA8l and IA9I) . Again, the (J x B) force cancels out 
since v \\ B. Combining cq l27l and [251 gives the momen- 
tum equation parallel to field lines: 



a 
V 



dv 
ds 



,d\i\B 

ds 



dU 
ds 



(29) 



At the critical point, v = a, and hence to avoid a diver- 
gent acceleration, the right hand side must go to zero, 
giving 



,d\nB dU 



(Is 



ds 



(30) 
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at the sonic point. The term on the left hand side rep- 
resents the pressure gradient due to the geometry set by 
the magnetic field. The term on the right hand side is 
the net acceleration along the field line. 

Since the sonic point is sufficiently close to the surface 
that the field geometry is not much perturbed by ex- 
ternal currents, we use dipole geometry to evaluate the 
sonic point position. The dipole field in cq[8] gives the 
intermediate result 



d In B 3 cos 8 



ds 



Nr 



1 



sin 2 
SN 



(31) 



The two terms on the right hand side of ea l31l repre- 
sent the field line divergence due to the r -3 factor from 
dipole field geometry, and the ^-dependent factor TV. The 
term due to differentiating N is negligible for large loops, 
but becomes important for small loops near the equator. 
Plugging eq[TT] and [3T] into cql3"0l and eliminating 8 us- 
ing the field line geometry in eq |12[ we find the following 
equation to determine the sonic point r = r s : 



GAL 



Q 2 r- 



3a 2 

r 



(32) 



where TV = yl — 3r/4r oq . This equation, solely in terms 
of r. can be solved as a function of the parameters A, e 
and r eq /R = 1/ sin 2 0&. 

First we examine the limit in which tidal forces can be 
neglected. This simple case highlights the importance of 
the magnetic field geometry, and would apply for slow 
rotating planets distant from the star. Setting f2 = in 
eq 1321 the simpler equation 



GM P 
3a 2 r 



= 1 



8r cqV /f - 3r/4r, 



(33) 



cq 



results. For large field lines r cq 3> R, field line curvature 
is negligible and the sonic point sits at 



GM P 
3a 2 



R 



(34) 



which differs from the spherical wind result by the 3, 
instead of 2, in the denominator. Including finite sin 2 8b, 
the sonic point moves inward somewhat. 

Next, we include tidal forces, but ignore the field line 
curvature terms on the right hand side scaling as r~^, 
a good approximation for field lines near the pole. This 
approximation eliminates the possibility that the tidal 
force can point outward. In this case, the sonic point 
equation becomes 



GAL 



3a 
r 



(35) 



The key point is that now the effective gravity on the left 
hand side has a minimum. If the pressure term on the 
right hand side is smaller than this minimum, a transonic 
solution is not possible. The solution of the cubic eg 1551 
lying near the planet disappears for sufficiently strong 
tidal forces 



i2 > f2 cr it — 2 



2tt 



GAL 



p 



2 a 

3 r sQ 



3.4 days 



0.7 M 



10 km s" 



(36) 



0.1 - 



10" 




FlG. 6. — Values of e versus \ for the transiting planets. The 
points are the data, with symbol type as in Figures [5] and [3] The 
line is the critical tidal strength e = 4/ A 2 above which the wind is 
suppressed in the polar region. We have used a = 9.3 km s _1 to 
make the plot. 

Eq |36l can also be written in dimensionless form as 



4 
A 2 



(37) 



When eq[3S] is satisfied, for planets sufficiently close to 
the star, there is no sonic point solution, and the wind is 
shut off near the poles. The density distribution on these 
field lines will be hydrostatic. Hence, for sufficiently 
strong tides such that the rotation velocity ri C rit^sO a t the 
sonic point is supersonic, a second dead zone is created in 
the polar regions where tides point inward. Whereas the 
first dead zone in the equatorial region is due to magnetic 
pressure dominating gas pressure, the second dead zone 
near the poles arises due to the large potential barrier. 

Figure |5] shows e versus A for the observed transiting 
planets. Except for a handful of planets with the smallest 
values of both e and A, most of the planets are in the 
strong tide limit with e > e cr it- The planets in the upper 
right hand corner will have the polar wind partially shut 
off, while the planets in the lower left hand corner will 
be able to drive a polar wind. 

Next, we retain the terms due to the tidal force in 
cq|3"2l Taking the limit Cl — > oo in eq|3"2l the sonic point 
in the strong tide limit is 



2r, 



■if 



(38) 



In this limit, the sonic point occurs at a fixed fraction 
of the loop equatorial radius, dependent only on cos 2 (f>, 
and agrees with the magnetic Roche lobe radius found 
in cq lI4l where the net force first points outward. 

Figure [7] shows an example numerical solution of eq |32l 
for the sonic point radius as a function of footpoint angle 
&b- Dipole geometry was used to produce this plot. In 
the weak tide limit (e — > 0), the solutions asymptote 
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Fig. 7. — Sonic point radius as a function of footpoint co-latitude 
for A = 12, P = 1 and a range of e (dashed lines). The three solid 
lines label the e = sonic point, r a o/R = A/3, the e = 4/A 2 = 0.16 
sonic point, r s /R = A/2, and the e = oo sonic point, r s = 2r eq /3/. 
Each dashed line terminates at large 9 at the dead zone, 9 = 9&. 
The dead zone shrinks (9^ 7r/2) as e increases. The value of /3 is 
needed only for the size of the dead zone. The longitude cos 2 <j> = 1, 
along the star-planet line, has been assumed here. 



to eq[M] for large loop size, and decrease slightly before 
terminating at the dead zone, 8 = 9d- For small e < 
4/A 2 , the sonic point moves out in the polar regions, and 
inward closer to the equator; the dividing line between 
these two behaviors depends on if the net force is outward 
or inward near the looptop. Next, for the critical value 
e = 4/A 2 , the sonic point is at roughly r s ~ (X/2)R near 
the pole. For larger values e > 4/A 2 , the sonic point 
near the pole jumps out to a radius much further from 
the planet, near r s ~ 2r eq /3/, where the net gravity 
changes sign. For small values of 9b, this solution may 
be at tens to hundreds of planetary radii, and is of no 
physical interest. Physically, when the sonic point moves 
outside the region of interest the field line is effectively 
hydrostatic. 

Figure E] shows the velocity at the base, Vb, for the 
same parameters as in Figure [7J This velocity was found 
by using the Bernoulli constant evaluated at the sonic 
point and the base. For weak tides (e = 0), the velocity 
at the base is nearly constant over the wind zone. As 
e increases, the base velocity decreases near the pole, 
where the sonic point has moved outward, and increases 
closer to the equator, where the sonic point has moved 
inward. 

At sufficiently large radii, pressure gradients rapidly 
become negligible and the fluid should move on a nearly 
ballis tic t rajectory. Using the Bernoulli integral defined 
in eq |A9[ with the sonic point as reference location, the 
velocity parallel to the magnetic field is 



1 2a 2 In ( ) +2(U S -U) 



B a 



rametcrs as Figure [7] 

where the subscript "s" refers to the sonic point. The 
transit signal depends only on the gas velocity within 
the stellar disk at r < R+, where i?* is the stellar radius. 
The tidal potential term — 2U dominates at large radius. 
For the purposes of a simple estimate, this asymptotic 
expression, evaluated at the stellar radius, gives 



ap ~£lR* (/sin 2 6-1) 1/2 
( 3.6 day 



: 24 km s" 



Uo 



/ sin 2 9 - 1 



1/2 

(40) 



(39) 



We caution the reader that the (logarithmic) enthalpy 
term in cq |39l is not negligible, and can increase the 
velocity at r = i?* by a factor of order 2 (sec Figure 
ITT]) . Though supersonic, the velocity in eq|TD] is much 
smaller than the ~ ±100 km s~ x at which absorption 
is observed in the Lyman a spectrum of HD 209458b 
(jVidal-Madjar et al.ll2003l ). and hence velocity gradients 
cannot be the origin of the observed transit depth. 

In the following sections, we will implement the general 
theory that describes the structure of the dead zone (§[6]) 
and of the wind zone (§ [7]) to construct global models of 
hot Jupiter magnctospheres for comparison with transit 
observations of HD 209458b. A simplified ID thermal 
model motivates our choice of the pressure and sound 
speed at the base of the global models, which are key pa- 
rameters for determining the magnctospheric structure. 

8. THE H AND H+ LAYERS: A SIMPLIFIED ID THERMAL 

MODEL 

The thickness of the warm H layer with tempera- 
ture T ~ 5,000 — 10,000 K is a crucial parameter in 
determining the transit de pth, as emphas i zed in the 
phenomenolog i cal m odel of iKoskinen et all (|2010l) (see 
I Garcia Mu noz ( 2007j) for a discussion of the role of base 
pressure for an atmosphere undergoing energy-limited es- 
cape) . The large temperature and small mean molecular 
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weight increase the scale height to ~ 0.1i? p h- If the warm 
layer extends over > 10 scale heights, the transit radius 
can be significantly increased over the photospheric ra- 
dius of the optical continuum. In this section, we con- 
struct a simple ID model in photoionization and thermal 
equilibrium to determine the depth of the warm H layer. 

At sufficiently low density and temperature, the rates 
of collisional ionization and 3-body recombination are 
slow compared to photoionization and radiative recom- 
bination, respectively, and the ionization state is set by 
a balance of the latter two processes. We consider a 
pure hydrogen gas for simplicity. Let n e , n p and rift be 
the density of electrons, protons and hydrogen atoms. 
Charge neutrality implies n e = n p . The ionization frac- 
tion in the "on the spot" a pproximation is found by soly - 
ing the algebraic equation (jOsterbrock fe Ferlandll2006|) 



n H J(N H ) =a B (T)n e n p = a B {T)n p . 
At the 50% ionization point, 

J 

;q — • 
OtB 



nn = n 



p 



(41) 



(42) 



Here a B (T) ~ 2.6 x 10~ 13 cm 3 s^lO 4 K/T) 8 
is the case B radiat ive recombination rate 
(jOsterbrock fc Ferlandl[2006l ). 



J(N H ) = / dv <p„ e- NH ^a pi (u) 



(43) 



is the photoionization rate per H atom, ip v is the pho- 
ton flux per unit frequency interval, Nft is the atomic 
hydrogen column from the point in question to the star, 
cr v i(v) = cr p i(i/ /^) 3 is the H atom bound-free cross sec- 
tion, the threshold cross section is er p i = 6.3 x 10 -18 cm 2 , 
and the threshold frequency is vq = 13.6 cV/h. The ex- 
ponential factor in eq |43l takes into account attenuation 
of the stellar radiation. Eqf43] may be computed as a 
functi on of Nh, as described in lOsterbrock fc Ferlandl 
(120061) . 

For the thermal balance, the dominant processes 
are photoelectr ic heating from the ionization of hy- 
drogen atoms (jYellel I2004D . and cooling by collision- 
ally excited Lyman a em ission from electron impacts 
(|Murrav-Clav et al.ll2009| ). Assuming 100% efficiency of 
turning photoelectron energy into heat, the heating rate 
per reaction is 

/>oo 

Q(N H ) = dv Vv e-tiB'A") api {v) h(v - Vq )(U) 



Balancing photoelectric heating and cooling by collision- 
ally excited Lyman a emission implies 



n H Q(N H ) = A(T)n e n H . 

Note that Uh cancels out of cql45l 
cooling coefficient for Lyman a is 
2.9 x IP" 19 erg cm 3 s " 1 



The 
A(T) 



(45) 
line 



_^10 4 K/Texp(-118,400 K/T) 
(Dalgarno fc McCravl 119721 ). Increasing distance from 
the star and decreased heating efficiency act to decrease 
Q. 

Eq. rjj] an d US] are two algebraic equations which can 
be solved for rift and n p in terms of Nh and P = 
(2n p +nH)kbT. In practice, we assume a trial T, compute 
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Fig. 9. — Photoionization rate (top) panel and heating rate 
(bottom) as a function of neutral hydrogen column, for a planet at 
D = 0.05 AU around a solar-type star.. 

rift and n p from ea l411 and then compute the imbalance 
of heating and cooling in cq |45l The temperature is iter- 
ated until thermal balance is achieved. The equations for 
dependent variables Njj and P and independent variable 
r are the definition of column 



dN 



H 



dr 



rift 



and hydrostatic balance 



dP 
dr 



GM p m p 
~2 ( n H 



(46) 



(47) 



where tides have been ignored for simplicity. The inner 
boundary condition is r = R at the chosen base pressure. 
The column Njj should go to zero at the outer boundary. 
We enforce this boundary condition at the finite, but 
large, radius r = 40i?. 

To compute the integrals in eqf43] and cqf44l we 
use the quie t sola r Lyman continuum spectrum from 
IWoods et al.l (|1998T ) , which tabulates J dvLp v in each fre- 
quency bin. The results are shown as the solid lines in 
Figure O Approximate power-law fits are also shown, 
along with curves representing pure exponential attenu- 
ation for comparison. A similar fit with shallower slope is 
shown in the bottom panel for the heating rate Q(Nft). 
At small optical depth, Jo ~ 6 hr and the mean pho- 
toelectron energy is Qo/Jq — 2.7 eV. Given the form of 
the integrand in eq |43l one might have expected an expo- 
nential scaling of the form J oc exp[— (constant) CpiNn], 
implying negligible heating and ionization deep in the H 
layer. This is shown as the dot-dashed line in Figure |9l 
and cuts off much too sharply. The numerical result is 
better fit with a power-law, J cx N^ 1 ' 5 , leading to larger 
heating rate deep into the H layer. 

The weak scaling of J with Njj can be explained as the 
competition between exp[— NH<Tpi{v)], which increases to 
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be solved analytically as 



red = 0.1AU 
black = 0.05AU 
blue = 0.025AU -3 




10" 3 10" a 0.1 1 

P[dyne cm -2 ] 

Fig. 10. — Temperature (top panel), number densities (second 
panel), neutral hydrogen column (third panel), and radius (bottom 
panel) versus pressure for a M p = O.TMj, R = 1.3Rj planet. The 
red, black and blue lines are for planets at D = 0.1, 0.05, 0.025 AU, 
respectively. The hollow circles show the position where Lyman 
continuum radiation (at threshold) from the star becomes opti- 
cally thick. The crosses show the 50% ionization point. The filled 
triangles show the position in the atmosphere where stellar Lyman 
a photons at ±100 km s from line center become optically thick. 



higher frequencies, and (p v , which on average decreases 
to higher frequency. Ignoring lines in y>„, the product of 
these two functions is a sharp peak, similar to the Gamow 
peak in thermonuclear reaction rates. For instance, if 
one approximates the Lyman c ontinuum as a blackbod y 
with temperature T = 8300 K (jNoves fc Kalkofenlll970l ). 
steepest descent evaluation of eq03] gives the weak ex- 
ponential scaling J cx exp[— 15.9(A r ^fcr p i) 1 / 4 ], which in- 

1 /4 

volves Nfj in the exponent rather than Nh- However, 
the blackbody fit is not adequate, as even this weaker 
scaling cuts off too fast. We find a better fit is to use 

a power-law form ip^ oc z^~ 7 , leading to J cx N H ^ +2 ^ 3 , 
with no exponential scaling. Choosing 7 = 2.5 then re- 
covers the observed scaling. While these analytic scalings 
are useful for intuition, the solar Lyman continuum con- 
tains many strong lines, and is not well approximated by 
a smooth continuum function when computing J. 

Figure [10] shows a numerical integration of eqs. I4"6]l47l 
I41l and l45l Spherically symmetric irradiation is assumed, 
the base pressure is chosen to be P = 1 dyne cm~ 2 , and 
the Nh — > condition is applied at the (arbitrary) ra- 
dius r = 40-Rph- Parameters appropriate for HD 209458b 
have been chosen with M p = 0.7Mj and R = 1.3Rj. In 
each panel, the three different color curves arc for differ- 
ent semi-major axis; the semi-major axis of HD 209458b 
is D ~ 0.05 AU. In the top panel, note that T increases 
slowly as the planet is moved closer to the star. A factor 
16 increase in stellar flux translated into only a 10-20% 
increase in the pressure range of interest, due to the ex- 
ponential T dependence in the cooling rate. Eq|45] can 



1.1 x 10 4 K 



1 + 0.089 ln[(n e /n eq , )(Qo/<2)]' 



(48) 



This simple estimate shows that the temperature is al- 
ways near 10 4 K near the base of the photoionized layer 
where nn = n p ~ n cq . As n e decreases below n eq , the 
temperature rises logarithmically slowly. This temper- 
ature rise toward small density will increase the scale 
height and density at larger radii relative to an isother- 
mal model with the same base conditions. 

The second panel shows the rapid inward increase of 
nn, with a change in slope at the H-H + boundary. The 
slow decrease of n p into the atmosphere is due to the slow 
decrease of J with Nh- If we had used the exponential 
scaling J oc exp(— NnCpi) shown in Figure [HI T and n p 
would have decreased inward much more rapidly. The 



hollow circles show the position of nn 



at each D. 



Planets further from the star remain neutral higher up 
in the atmosphere. 

The third panel shows Nh, and x's show the position 
of Nh<J v \ = 1, where the Lyman continuum at threshold 
is optically thick. The optical depth unity point moves 
to lower pressure for planets more distant from the star, 
although the radius at this point is relatively constant. 

The bottom panel shows radius in units of base radius. 
The warm H layer with temperature T ~ 10,000 K at 
pressures P = (0.1 — 100) nbar contributes an amount 
~ (1 — 2) x R p h to the radius. Specifically, even the region 
below Nn<Jpi = 1 can be warm enough to contribute 
significantly to the radius. In our reference global model 
below, we use a base pressure P ss = 40 nbar (see Model 
1 of Table 1). 

We now discuss how the simple ID model differs from 
previous investigations. One crucial difference is that 
the dead zone should be hotter than the wind zone, 
for which a diabatic expansion is an important coolant. 
lYelld ()2004h included heating arising from photoioniza- 
tion of hydrogen, but ignored the attenuation of the stel- 
lar Lyman continuum into the atmosphere, leading to 
an overestimate of the heating rate. This attenuation 
of EUV in the H layer will also lead to smaller heating 
by H2 photoionizat i on m uch deeper in the atmosphere. 
iMurrav-Clav et all ()2009[ ) included finite optical depth, 
but enforced an exponential cutoff which led to a steep 
temperature drop below the NnVpi = 1 point. As shown 
in Figure [5] the exponential cutoff is too rapid, and the 
slower power-law cutoff found here gives additional heat- 
ing deeper in the atmosphere. 

Figure [10] shows that the tin = n p and Nn&pi = 1 
points occur near each other at D = 0.05 AU, hence 
attenuation due to finite optical depth cannot be ignored. 
The relative position of the tih = n p and Nh<J p1 = 1 
layers can be estimated as follows. The number density 
at which n v = tih is 



J 

( T 
1.8 x 10 8 cm -3 ' 



VlO 4 K 



'0.05 AU X 2 



D 



9.) 



This can be converted into a pressure as 
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8.0 x 1CT 4 dyne cm" 2 

2 



Mj 



10 



id 



R 



0.05 au 
D 



T 



10 4 K 



0.8 



(50) 



The D~ 2 scaling implies that the atmosphere becomes 
neutral out to smaller pressures as the planet moves 
away from the star. To estimate where Ng o~ v \ = 1, 
we assume n p 3> riH, P — 2kbTn p , and eq|4ll gives 
TiR — n 2 /n cq ,o, where n cqj o = Jo/as- For H atom scale 
height ~ kbT/m p g, the pressure at which Nh&pi = 1 is 



Pi 



:2k h T 



kbTapi J 



1/2 



oc D~ 



(51) 



The scaling with D in ea l51l implies that planets further 
from the star will become optically thick at lower pres- 
sure. Comparing the scalings in ea l50l and 1511 one may 
expect 50% ionization to occur outside Nh<J p1 = 1 at 
sufficiently large D. 

Lastly, we note that since the cooling rate oc n e riH, 
we may expect Lyman a cooling to be inefficient at both 
high and low density, where other cooling mechanisms, 
such as heat conduction, may be more efficient. 

9. GLOBAL MODELS 

We now attempt to construct global models for the 
magnetosphere, in order to compute the planetary trans- 
mission spectrum and mass loss rates. We assume the 
follow ing magnetic field model (|Mestell 119681 : lOkamotol 
IT971 : 



S (|) (e r cos0+ |e e sin0) , (r < r& ) 
B{r,0)={ Kr ' ( n3 2 ; [52 

Bo(£) (^) 2 cos0e r , (r>r d ). 

This global field model allows a position (r, 9) to be as- 
sociated with a base co-latitude 9b at r = R. The field 
is dipole for r < r^, and radial outside r^. This radial 
field is distinct from the split monopole due to the cos 9 
factor. Eq J52I is approximately correct for the dead zone, 
and also for determining the sonic point position in the 
wind zone if the sonic point is close to the planet. How- 
ever, as ballistic trajectories with speeds far less than 
escape speed are expected to be bent down toward the 
orbital plane, the radial field assumption will be unphys- 
ical for some latitudes. 

Given parameters A, e, (3, and longitude tp, we first 
solve eq |26l for rd using dipole geometry. For field lines in- 
side the dead zone, sin 9b > sin 9d = yRjvd-, the velocity 
is zero and the density is given by the hydrostatic expres- 
sion in eq |19t where P ss is a model parameter. On field 
lines outside the dead zone, sin 9b < sin 9d = y/R/ra, we 
search for sonic points by looking for minima of the quan- 
tity — ln(£>/_B ) — Ufa? on field lines, between the base 
radius R and a chosen maximum radius r max . Given the 
position of the sonic point (r s , 9 S ), the Bernoulli equation 
v 2 /2 + a 2 \n(B/v) + U = constant can be used to solve 
for the velocity at the base, Vb(9b,4>). Given Vb and the 
base density (eq ITS)) , the Bernoulli equation may again 
to used to find the run of v and p on the field line. 

Detailed results will be presented for the 9 models 
listed in TablcQ] The planetary mass and radius, and the 



stellar mass and radius are characteristic of HD 209458b. 
We vary parameters not directly measured, such as P ss , 
a, Bq, as well as the orbital radius D. 

The numerical implementation of the sonic point solver 
deserves further discussion. Anywhere from zero to sev- 
eral solutions to eg 1301 may be found in the interval 
R < r < r max . Some sonic points may be spurious if 
a potential barrier exterior to the sonic point decelerates 
the flow to subsonic, even zero, speed. These spurious 
solutions are discarded by defining the true sonic point 
solution to be a global minimum of — ln(_B/i?o) — U/a 2 
which occurs within the interval R < r < r max . If the 
global minimum occurs at either of the endpoints of the 
interval, then there is no good sonic point solution. 

For minimum at the base r = R, the integration is 
flagged as a region of parameter space with no solution, 
as we should have used a base position deeper in the 
planet. For sonic points sufficiently deep in the planet 
that ram pressure dominates magnetic pressure at the 
sonic point, we expect the Roche lobe overflow model to 
be recovered. The other problem is that the global min- 
imum can occur at the outer boundary of the integra- 
tion, r = r max . For instance, this can occur in the polar 
regions due to the upwardly increasing potential. For 
r > rd, field lines with sin# > sin# cr i t = / _1 ^ 2 have out- 
ward tidal force. Field lines with outward tidal force will 
have fluid accelerated outward, promoting the existence 
of a sonic point. The field line starting at base co-latitude 
sin9b,crit = (R/ ffd) 1 ^ 2 will be the last field line on which 
the tidal force points outward. Accordingly, if no sonic 
point solution is found in the radial range R < r < r max , 
and the field line has sin# < sin 9b. cr it, then we treat 
the field line as hydrostatic and set the velocity to zero. 
Such field lines have outwardly increasing potential, and 
hence outwardly decreasing density. If, however, we had 
used a field line model that allowed the field at r > rd 
to bend downward toward the orbital plane, it is likely 
that a sonic point could have been found. This affects 
our later numerical results, as we analytically predicted 
the critical tidal strength to be e ~ 4/ A 2 in order to find 
sonic points in the polar region, whereas the above pre- 
scription would force these field lines to be hydrostatic 
due to the potential barrier. This approximate treatment 
of the polar regions likely does not affect either the total 
mass loss rate, or the column density profiles, since the 
sonic point will be so far from the planet that the veloc- 
ity near the planet is quite subsonic, and the fluid will 
be nearly hydrostatic. 

Figure [TTJ shows the density and speed along a partic- 
ular field line in the wind zone. The change in slope near 
r = 4.7 R is the change in field geometry at r — rd- The 
speed along the field line is somewhat larger than the 
asymptotic result in cq 1401 due to the enthalpy (In) terms 
in eq!391 In the lower panel, the density is approximately 
hydrostatic inside the sonic point at r s = 2.4i?, and de- 
creases roughly as p oc B/v oc r~ 3 outside that point. 
This plot explicitly demonstrates that the gas density in 
the wind zone can be orders of magnitude smaller than 
nearby gas in the dead zone, which satisfies hydrostatic 
balance. 

Figure [T^] shows contours of mass density on slices 
through the center of the planet in the y — z plane, as 
viewed during transit, and the x—z plane, as viewed mid- 
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TABLE 1 

Global Models 



model 




R , b 


D c 


R d 


1 ss 


a f 




A h 


i 


M 


' L 


l 


' a 


SF/F n 


M° 


' scat 


# 


(Mj) 


(Rj) 


(AU) 


(Rj) 


(^tbar) 


(km/s) 


(G) 








(R) 


(R) 


(R) 




(g/s) 


(R) 


i 


0.7 


1.3 


0.05 


1.4 


0.04 


10.0 


8.6 


9.3 


0.032 


0.054 


4.6 


3.1 


6.0 


0.074 


9.3 x 10 10 


6.6 


2 


0.7 


1.3 


0.05 


1.4 


0.4 


10.0 


8.6 


9.3 


0.032 


0.54 


4.6 


3.1 


3.3 


0.31 


1.4 x 10 12 


10.0 


3 


0.7 


1.3 


0.05 


1.4 


0.004 


10.0 


8.6 


9.3 


0.032 


0.0054 


4.6 


3.1 


9.7 


0.025 


6.6 x 10 9 


4.5 


1 


0.7 


1.3 


0.05 


1.1 


0.04 


13.0 


8.6 


5.5 


0.019 


0.054 


1.6 


1.8 


3.0 


0.24 


1.7 x 10 12 


9.7 


5 


0.7 


1.3 


0.05 


1.4 


0.04 


7.5 


8.6 


16.6 


0.056 


0.054 


1.6 


5.5 


23.2 


0.030 


1.5 x 10 s 


2.8 


6 


0.7 


1.3 


0.05 


l.i 


0.04 


10.0 


43.0 


9.3 


0.032 


0.0022 


1.6 


3.1 


11.6 


0.091 


5.9 x 10 10 


6.8 


7 


0.7 


1.3 


0.05 


1.4 


0.04 


10.0 


2.9 


9.3 


0.032 


0.49 


1.6 


3.1 


3.4 


0.072 


1.4 x 10 11 


7.8 


8 


0.7 


1.3 


0.025 


1.1 


0.04 


10.0 


8.6 


9.3 


0.25 


0.054 


2.3 


3.1 


6.1 


0.075 


6.6 x 10 1() 


1.2 


9 


0.7 


1.3 


0.10 


1.1 


0.04 


10.0 


8.6 


9.3 


0.0040 


0.054 


9.2 


3.1 


5.9 


0.096 


1.4 x 10 11 


7.7 



a Planet mass. 

b Transit radius. 

c Orbital separation. 

d Radius at the base of the warm H layer. 

e Pressure at base radius (r = R) at the substcllar point (9, <f>) = (tt/2, 0). 
f Isothermal sound speed. 

g Magnetic field at pole (9 = 0,7r) at base radius (r = R). 
h A = GMp/Ra 2 

1 e = (QR/a) 2 , where fi is the orbital frequency. 

1 /3 = 8ir P BB / (Bq / 2) 2 is the base value at the substcllar point. 

k Radius of L\ — L2 Lagrange points. 

1 Sonic point radius ignoring tides and field line curvature (eq l34J l, 
m Dead zone radius in <j> = tt/2 plane. 

n Integrated Lyman a transit depth from —200 to 200 km s _1 from line center. The disk inside the base radius is asssumed opaque, and 
contributes (R/R*) 2 ~ 0.015 to the transit depth for R = lARj and i?* = 1.15-R©. 
° Computed mass loss rate. 

p Radius (A Bca t /it) 1 / 2 corresponding to area A sca t in x — z plane over which t v > 1 for Lyman a cross section a v = 10~ 15 cm 2 . This cross 
section corresponds to a velocity ~ 25 km s — 1 from line center. 
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tity plotted is p/m p = nn+n p (note that this quantity is 
distinct from the total number density n to t = + 2n p , 
which depends on the details of the photoionization 
model). Model 1 parameters listed in Tabled] were used. 
Near the planet the contours are approximately spheri- 
cal, since the velocities are everywhere subsonic and the 
tidal force is small. The bulge at the equator is the equa- 
torial dead zone. The poles are hydrostatic as tides have 
shut down the wind, and the inward tidal force at the 
pole causes the density to decrease outward faster than 
at the equator. The impact of the tidal force on the dead 
zone can be seen by comparing the upper and lower pan- 
els in Figure [121 Along the x-direction, the outward tidal 
force decreases the size of the dead zone, but the same 
tidal force also causes the dead zone to have higher den- 
sity. Outside the sonic point, the density in the wind 
zone is smaller than in the neighboring dead zone, which 
pushes the density contours inward in the wind zone. 
Near the equator, outside the dead zone, the density be- 
comes quite small, hence the pile-up of contours near the 
critical angle sin0 cr i t . 

In the next section, the wind models in Table 1 are used 
to compute the planetary mass loss rate in the wind zone. 



Fig. 11. — Velocity and density versus radius along a field line 
in the wind zone. Model 1 parameters from Table [T] are used. The 
field line is located at cf> = with 9 b = 0.35 rad. The dead zone is 
at 9 d = 0.48 rad and the critical field line is at t9 b crit = 0.23. The 
discontinuity at r/R = 4.7 is the dead zone radius, where the field 
changes shape. T he l ine labeled "tide only" is the asymptotic ap- 
proximation in co 1391 and the line labeled "hydrostatic" evalu ates 
the density using the hydrostatic balance approximation in eq |19l 



way between primary and secondary transit. The quan- 



10. MASS LOSS RATE AND SPIN-DOWN TORQUE 

The mass loss rate is computed by integrating pv r = 
pb r vi, over the surface area of the wind zone at the base. 
Using the base density from cq[l8l the mass loss rate is 



M = R 2 Ps 



ddb sin 6>b b r Vb 



wind 



Pb{0b,<j>) 
Pss 
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Fig. 12. — Contours of mass density, log 10 p/m p [cm -3 ] in the 
y — z plane (upper) at x = 0, and in the x — z plane (lower) for 
y = 0. Model 1 parameters are used (see Table [TJ. 

= 4.0xl0" g3 - ( T A-) 3 
where the dimensionless integral 

.F(W)=8 f% [ 9d d9 b sin6 b U*) (^4^) 
Jo Jo \aJ\ p ss J 

The mass loss rate for fixed M, R and Bq, but varying 
a, D and P ss is shown in Figure [T3] The steep decline 
of M with A is due to smaller density at the sonic point 
radius. The mass loss decreases slightly for large e due 
to the smaller fraction of open field lines. 

Why is the mass loss rate M proportional to the base 
pressure P ss ? Recall that we are approximating the true 
atmosphere with an isothermal model. The appropriate 
values of a 2 and P ss , as determined by photoionization 
equilibrium and heating/cooling balance, have been dis- 
cussed in section[5J The sonic point lies at a fixed radius, 
given roughly by eq|M]bascd on the choice of sound speed 
a. The location of the base of the isothermal layer is also 
at a fixed radius, estimated to be l.li? p h (see §[4}. By 




Fig. 13. — Total mass loss rate as a function of A for M p = 
0.7Mj, R = lARj and B = 8.6 G. The upper (lower) set of three 
curves uses base pressure P ss = 0.04 (0.004)/xbar. The line style in 
each set of three curves gives the value of e, the tidal strength. 

cq[T9l p(r = r s ) cx p 8S , therefore the density at the sonic 
point is proportional to the base density p ss (and there- 
fore also P ss )- Consequently, if a larger value of P ss is 
required to explain the transit depth, the mass loss must 
be increased proportionally. 

The mass loss rates in Figu re [T5| are largely consis- 
tent with previous studies (e.g.. lMurrav-Clav et al.ll2009f) 
when comparable gas density is used. By compari- 
son, an unmagnetized, spherically symmetric, is other- 
mal wind would have (jLamers fc Cassinellil [19991 T ~ 
7rA 2 exp(3/2 — A), which would be a factor of ~ 3 — 10 
larger than the curves in Figure I13[ and with a slightly 
flatter slope. Inclusion of the magnetic field decreases 
the mass loss rate, mainly due to the decrease in area 
occupied by the wind zone. 

The angular momentum loss rate depends on the ra- 
dius at which the torque is applied. For an isolated 
planet, the field lines remain rigid out to the Alfven ra- 
dius. But this location may be at many tens of planetary 
radii, and may be pre-empted by the interaction of the 
planetary wind with the stellar wind. By assuming the 
torque is exerted at a radius rtorquc we estimate an an- 
gular momentum loss rate 

A^r t 2 orquo -7.3xl0 28 erg 

x (i^Si) (lo^gs- 1 ) ' (55) 

While this torque may cause moderate changes in the 
spin rate for an isolated planet on Gyr timcscales, it 
likely not large enough to torque the planet away from 
synchronous rotation to the extent tha t significant grav- 
itatio nal tidal heating will occur (see lArras fc Socratei 
120091 for a discussion of the necessary torques). 
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11. NEUTRAL HYDROGEN COLUMN DENSITY 

Given the global MHD models derived in section^ we 
require a model for the run of ionization in the magneto- 
sphere in order to compute observable quantities such as 
the transmission spectrum. In section [5J we discussed a 
model including only the dominant processes: photoion- 
ization and radiative recombination of hydrogen. In the 
remainder of the paper, we use the simpler optically thin 
limit to evaluate nn given p: 



'cq,0 



V^cq.O 



(56) 



Here n eq ,o = Jo/&b is the density at which Hh = n p in 
the optically thin limit. The use of Jo instead of J simpli- 
fies the calculation, as only the local density is required, 
and not the column Nh- This approximation may un- 
derestimate riH near the H-H + transition, but should be 
adequate for our purposes. 

To evaluate the neutral hydrogen column density, 
we first evaluate p (§ EJ) , and then njj (eqJSBJ), on 
a grid of [x, y,z), with each coordinate in the range 
(— l.LR*, l.liJ*). The column is displayed as seen at 
transit, i.e. we integrate over the coordinate along the 
star planet line to get column 



N H (v,z) 



l.ifl. 



dx nn{x, y, z) 



(57) 



-1.1R* 



as a function of the impact parameters y and z. Figure 
1141 shows contours of hydrogen column density for the 9 
models listed in Table [TJ All models have M p = 0.7 Mj 
and 7? = lARj, and vary a single parameter P ss , a, 
Bq and D in turn. The fiducial case, Model 1, clearly 
shows the equatorial and polar dead zones, as well as the 
mid-latitude wind zone with comparatively smaller Njj. 
Model 2 (3) has P ss larger (smaller) by a factor of 10. 
This has the effect of decreasing (increasing) the dead 
zone size as well as scaling up (down) the density in the 
dead zone. Model 4 (5) has larger (smaller) a, leading 
to larger (smaller) density at a given distance from the 
planet, as well as increasing (decreasing) the size of the 
dead zone. Model 6 (7) has larger (smaller) B , which 
increases (decreases) the size of the dead zone. Model 8 
(9) has smaller (larger) D. Larger tide is more effective 
in shutting down the wind at the pole, but also decreases 
the size of the dead zone. 

Aside from the overall magnitude mainly set by the 
base pressure P ss and sound speed a, the dominant pa- 
rameters determining the appearance of each plot arc 
the equatorial and polar dead zone sizes. The equatorial 
dead zone size (see Figure [5]) depends on A, /3 and e. The 
size of the polar dead zone is set by the strength of the 
tidal force. The critical tidal strength in eq[37] refers to 
shutting down the wind at 9b = 0, and assumes dipolc 
field geometry. For e > e cr it, a range of 0b near the pole 
can have the wind shut off. The result depends on the 
which field geometry is chosen. For instance, the sonic 
point eg 1321 can be rederived for radial field lines. The 
discriminant of this cubic equation can be used to show 
that no sonic point can be found for 
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The critical tidal strength for radial field lines is e cr ;t = 
32/(27A 2 ), a slightly different numerical coefficient than 
the dipole case. As e increases above e C rit, the size of the 
polar dead zone increases. In the limit e 3> e c riti the wind 
is shut down in the entire region sin 2 6 < 1/f where the 
tidal force is inward. For large e and small /3, the polar 
and equatorial dead zones can dominate the volume near 
the planet (e.g., Model 6). 

12. LYMAN a TRANSMISSION SPECTRA 

In section [TT] we focused on understanding the hydro- 
gen column as a function of impact parameter, including 
the dependence on unknown parameters such as temper- 
ature and magnetic field. An additional effect on the 
transmission spectrum is the velocity gradients in the 
wind, which were studied in sections [7] and [HI In this 
section we compute the Lyman a transmission spectra 
for the global models, including both column and veloc- 
ity gradient effects. 

The transmission function, T Vl is the fraction of stellar 
flux at frequency v which passes through the planet's 
atmosphere without suffering scattering out of the beam. 

In terms of the out-of-transit stellar flux, F„°\ and the 
in-transit flux, F Vl T v is defined as 



T v = ■ 



1 V 



(0) 



(59) 



If the interstellar medium (ISM) optical depth, Tj) ^ , is 
constant over the stellar disk, and in time, and the geo- 
coronal emission is independent of time, then the ratio in 
eg 1591 depends solely on the properties of the planetary 
atmosphere, and is the fundamental quantity to compare 
to the data. 

We compute T u as follows. Let <j v be the Lyman a line 
cross section. We simplify the problem by assuming the 
planet to be at the center of the stellar disk. The optical 
depth through the planet's atmosphere at position (y, z) 
on the stellar disk is 



r (p), 



(y,z) = / dx n H {x,y 1 z) a v (x,y,z). 



(60) 

Assuming the stellar intensity is uniform over the disk, 
1 



dydz e 



>{y,z) 



(61) 



where the integral extends over y 2 + z 2 < i? 2 . As 
integrated measure of the transit depth, we compute 



~F 



(62) 



where ij)*^ and r^ ISM ^ are the unabsorbed stellar in- 
tensity and ISM optical depth, both assumed uniform 
over the disk. In practice, we follow IBen-Jaffell (|2008h 
and integrate over —200 km s _1 < Av < 200km s _1 . 
The interstellar medium (ISM) is assumed to have 
a temperature Ti sm = 8 000 K and hydrog en column 
AT H . isrn = 10 18 ' 4 cm" 2 ([Wood et all l200l . implying 



the line is dark inside lincwidth Av = c(y — vo)/vo < 

50 km s _1 . For 7^ we use the following (unnor- 
(58) malized) fit to the quiet solar Lyman a spectrum 
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Fig. 14. — Contours of hydrogen column density, log 1 Q(A r //[cm 2 ]), versus impact parameter in the y — z plane observed during transit. 
The parameters for each plot are matched to the labeled model numbers in Table [l] 



presented in | Feldman et al.l (I1997D (downloaded from 
|http : / / www.mps . m pg . de / pro j ects / soho/sumer /FILE/| 
Atlas.html): 



Av 



67 km s" 



(63) 



The Voigt function H{a,u) (jMihalas! [1971) is used for 
the line profile, giving 



o v = /12-t=t H(a D , u) 

m e C yJlTlXVD 



(64) 



where — e is the electron charge, m e is the electron mass, 
and f\2 = 0.42, A = 1215A and i> = c/A are the Ly- 
man a oscillator strength, line center wavelength and fre- 
quency. The Doppler width is Ai/p = VqV^/c where the 
hydrogen atom thermal velocity is v t h = (2fc;,T '/trip) 1 / 2 . 
The damping parameter is an = r/47rAz^D, where the 
natural linewidth is T = 6.25 x 10 8 s _1 . Finally, the dis- 
tance from line center, in Doppler widths, including both 



bulk motion and thermal broadening, is 

v - VQ Vx_ 

Av D w t h 



u - 



(65) 



where v x is the bulk motion directed from planet toward 
star, which is away from the observer. 

There are three instructive limits of cq(HI]to guide the 
intuition. First, if v x = and the gas is optically thick 
over an area A tran , with negligible optical thickness out- 
side this area, the fraction of flux absorbed by the planet 
is 



1-T„ = 



= 0.013 



tt(1.3Rj) 2 



1.15-R 



(66) 



Next, if v x — and the gas is optically thin, then the 
transit signal due to the optically thin area is 



i-r„ 



i 



dydzT<?\y,z) = {T<?)), (67) 



which is just the area-averaged optical depth, and is pro- 
portional to the total number of hydrogen atoms times 
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Fig. 15. — Lyman a cross section as a function of velocity width 
from line center. Both profiles are given for a temperature T = 10 4 
K. The solid line is for zero bulk velocity away from the observer, 
while the dashed line is for Ad = 50 km s . For bulk velocity 
toward the observer, the dashed curve would be reflected about 
Av = 0. 

their mean cross section. The third limit is when ther- 
mal motions are much smaller than bulk motions, and 
the line profile can be approximated as a delta func- 
tion S [v — Uq(1 — v x /c)]. In this case, the cross section 
is only nonzero at those values of = x*(v,y,z) where 
v x (x*,y, z) = c{vq — v)/vq is satisfied, so that the photon 
is shifted to line center in the atom's frame. In this case 
the optical depth becomes 



tI p) (v,z) = n H (x+,y, z)^— fu 



A 



:2xio _ 3 (n H (x*,yz) 
\ 1 cm d 



m e c \dv x (x ir ,y,z)/dx 
Porb \ n 



1 day/ \dv x (x*,y,z 



where in the second equality we have scaled the velocity 
gradient to the orbital frequency Q. Eq[68] shows that for 
hydrogen densities tih > 10 2 ~ 3 cm~ 3 , the optical depth 
along a line of sight will be high provided that there 
is gas with sufficiently large velocity to absorb at that 
wavelength. 

Figure [15] shows the cross section as a function of fre- 
quency in velocity units, at T = 10 4 K and for Aw = 
and Av = 50 km s — 1 . For HD 209458b, the transit radius 
is i? p h = 1.3Rj and the stellar radius is R+ = 1.15-R©, 
giving a transit depth SF/F = 0.013 in the optical con- 
tinuum. To explain the line-integrated Lyman a transit 
depth ~ 9% (e.g., see the discussion in iBen- Jaffelll2008() 
one could invoke an opaque disk of area ~ 7r(2.6i? p h) 2 . 
The central issue is that this disk must be opaque at 
Av > ±100 km s _1 from line center, requiring large 
columns of neutral hydrogen at radii 2 — 3i? p h- 

In Figure I10[ triangle symbols show where Lyman a 
radiation at frequencies ±100 km s" 1 from line center 
is optically thick on a radial line outward. This point 
is much deeper in the atmosphere from where Lyman 
continuum at threshold becomes optically thick, due to 



the rapid decrease in Lyman a cross section. Clearly 
in order to model the transit spectrum in the wave- 
length region of interest, one must include regions down 
to 1 — 10 nbar in the atmosphere. To quantify this 
statement, we compute the optical depth through the H 
layer where nu — p/m p is given by eq !17l Assuming the 
dominant contribution arises from the layer of steeply 
falling density, the slant optical depth is dominated by 
the region near x = and we find 



T<?\y,z)=a v 



p{x = 0,y,z) 



x / dxexp 



1 / GM p 
2^2 I 63 



3ft 2 x 2 



p(x = 0,y,z)\ ( 2iTb 3 /\R 



1 - {b/r L f 



1/2 



/100 km s- x \ 2 fP{x = 0,y,z)\ /10 km s" 1 
V Av J V 1 nbar 



1(A l/2 /6X 3/2 
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(1 - (Wr L ) 3 ) 



-1/2 



(69) 



Here b = \Jy 2 + z 2 is the impact parameter, cqj3] was 
used for the tidal potential, is given by eqfTl and the 
last equality assumes the cross section is on the damping 
wing (see Figure [T5]) . In the H + layer, eqU^l should be 
multiplied by 1/2 to account for the smaller H atom scale 
height. Eq[69] agrees roughly with the position of the 
triangles in Figure 1101 keeping in mind that the slant 
length is a factor of a few larger than the scale height. 
Eq l69l shows that the Lyman a transmission spectrum at 
Av = ±100 km s _1 is probing down to < nbar pressures, 
depending on the value of b/R. 

To give a more precise numerical estimate of the transit 
depth, we first compute the integrated quantity SF/F as 
in eq|55]for the 9 models in Table [TJ The result is given 
in the table. The velocity range is taken to be —200 < 
Aw[km s _1 ] < 200. Since 1 — T v decreases away from 

line center, (1 — T v )Fv°^ is peaked somewhat closer to line 

center than F„ , the amount depending on the details of 
the atmosphere. Transit depths of the correct magnitude 
SF/F ~ 5 — 10% can be achieved by adjusting the main 
parameters, P ss ~ 10 — 100 nbar and a ~ 8 — 12 km s _1 
to have values as expected from the ID model in Figure 
1101 The parameters Bo and D have a lesser impact by 
comparison. 

The frequency dependent transit depth, 1 — T v , was 
computed as in eqEH for the 9 models listed in Table 



[TJ and com pared to the data for (Fj — F v )/Fv ' from 
Figure 6 of iBen-Jaffell (|2008h . The results are shown in 
Figure 1161 Near line center, nearly the entire planetary 
atmosphere is optically thick, and absorption is nearly 
complete. Moving out from line center in the Dopplcr 
core, the cross section eventually becomes small enough 
that part of the atmosphere becomes optically thin, after 
which 1 — T v decreases rapidly. The curves level out 
when the damping wing is reached, after which 1 — T u 

decreases slowly as the = 1 point moves deeper into 
the atmosphere as P(x = 0, y, z) cx Av 2 . 

Given the large error bars, a range of parameter space 
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Fig. 16. — Frac tion al flux decrease, 1 — T„, versus frequency in velocity uni ts for HP 209458 b. Curves for Models 1-9 from Tablc[T]arc 
computed from cg l61l and points with error bars are the data from Figure 6 of Bcn-Jaffcl (2008). 



agrees with the data if the warm H layer extends suf- 
ficiently deep. For instance, Model 3 with base pres- 
sure P ss = 4 nbar is well below the d ata points with the 
smalle st error bars, in agreement with lMurrav-Clav et al.l 
(|2009l ). The most sensitive parameter dependencies are 
with the base pressure, P ss , and the sound speed (temper- 
ature), a. Increasing the magnetic field has the effect of 
increasing 1 — T v due to larger Nh in the magnetospherc. 
Somewhat offsetting this effect is that increasing B de- 
creases the size of the wind zone, which decreases absorp- 
tion near line center due to velocity gradients. Perhaps 
counter-intuitively, moving the planet further from the 
star increases the transit depth. Inspection of Figure [10] 
shows that the H extends to both lower pressure and 
larger radius for more distant planets with atmospheres 
in photoionization equilibrium. Lastly, we note that ve- 
locity gradients are only important for Av < 50 km s , 
and are more important for smaller D due to the larger 
tidal force. 

We end this section with a brief discussion of scatter- 



ing of Lyman a from H atoms in the magnetospherc. 
The problem with observing Lyman a during transit 

is that large Nh is required to create Tv ~ 1 at 
Av > 100 km s _1 . By contrast, at line center the cross 
section is ~ 10 5 times larger, implying the atmosphere 
is optically thick at line center out to much larger radii. 
We suggest that scattering of stellar Lyman a during the 
orbital phases in which the planet is moving toward or 
away from the observer may be detectable, and provides 
a probe of thermal gas in the magnetosphere, comple- 
mentary to the transmission spectrum measured during 
transit. During the orbit, the Doppler shift of the scat- 
tered spectrum varies in time due to the variation in 
line-of-sight orbital motion. The orbital velocities natu- 
rally produces a feature in the spectrum well outside the 
line core where ISM absorption dominates. 

For a planet in circular orbit, there is no relative radial 
motion with respect to the star, and the stellar spectrum 
at the planet is not Doppler shifted. However, when an 
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H atom in the planet resonantly scatters a stellar photon, 
that H atom is moving with respect to the observer due to 
the planet's orbital motion. Photons emitted by the star 
near line center (Av < 67 km s _1 ) have their frequen- 
cies shifted by w or b = 210 km s _1 (l day/Porb) 1 ^ 3 (for a 
solar mass star) due to the planet's orbital motion. To 
assess the area presented by the magnetosphere, we com- 
puted the area in the x — z plane for which r^ p ' > 1 for 
<r v = 10 -15 cm 2 , which corresponds to Av = 25 km s _1 
from line center for T = 10 4 K. The results are tabulated 
in Table [1] We find that the effective radius of the scat- 
tering disk is r sc ~ (5 — 10)i? for the models shown. The 
scattering disk for Lyman a is significantly larger than 
the radius inferred during transit. Assuming none of 
the resonantly scattered Lyman a photons are ab sorbed, 
and a lso assuming the Lambert phase function (jHapkei 
|1993| ) as an estimate, the reflected flux is F Te ^{v) = 
F*(i/')(2/3n)(r sc / D) 2 , where v ~ v ±v OT b is the observed 
frequency, v' ~ vq was the frequency emitted by the star 
before Doppler shift, and F^v) is the stellar Lyman a 
spectrum. The size of the reflected flux relative to the 
flux emitted by the star out on the wing at frequency v 
is then FreflOVi^O) = {F*{v')/F*{v))(2l^){r sc /D) 2 . 
Inspection of eq(HS] shows that the line center flux is ~ 
30(Aw/200 km s -1 ) 3 times larger than that at Av. This 
acts to enhance the scattered flux signal relative to the 
background flux level. Numerically we find the ratio of 
scattered, Doppler shifted flux to background stellar flux 
is then F refl (z/)/F*(z/) ~ 0.4(r sc /10Fj) 2 (l day/P orb ) 7 / 3 . 
While this signal may be small for HD 209458b at 
P orb = 3.5 days, F Ie& (v)/ F^v) ~ 0.02(r sc /10i?j) 2 , for 
planets with P al ±, = 1 — 2 days it may be large enough 
to be observable. 

13. COMPARISON TO ROCHE-LOBE OVERFLOW 

The magnetic wind model developed in this paper dif- 
fers in several res pects from purely hy drodynamic mass 
loss models (e.g.. iLubow fc Shulll975l ). In the standard 
Roche-lobe model for nearly equal mass stars, nearly all 
the gas leaves the donor in a narrow, cold stream through 
the LI Lagrange point. The first assumption underly- 
ing this solution is that r Sj o S> rt,i, so that the gas is 
subsonic at the LI equipotcntial for most (#,</>). From 
eq[7] and eq|34l this ratio is r^o/Vx = (eA 2 /9) 1 / 3 , and 
hydrodynamic Roche lobe overflow requires e 3> 9/ A 2 . 
Figure [6] plots e versus A, and shows that most, but 
not all, transiting planets are indeed in the r s ,o 3> tl 
regime; ignoring magnetic effects, Roche lobe overflow 
would then be a good approximation. In the opposite 
limit of e <C 9 /A 2 , the solution would more closely re- 
semble a thermally driven wind weakly perturbed by 
tides. The second assumption underlying a narrow flow 
through LI is that the mass ratio of the two bodies 
is near unity. Although the tidal expansion r <C D 
in eq[3] ignores the difference in potential between the 
LI and L2 Lagrange points, inclusion of higher order 
terms give s Ul2 — Uli — 2GMp/(3D ) for the potential 
difference (jMurrav fc Dermottl l2000t ). When the ratio 
2GM p /(3Da 2 ) = (2/3)(eAM p /M*) 1 / 3 < 1, the density 
difference between the LI and L2 points is small, and 
nearly equal mass loss is expected through LI and L2. 
While mass loss through LI enters into an orbit around 
the star, mass loss through L2 leads to gas in a circumbi- 



nary orbit. 

MHD effects, in particular the existence of a dead zone, 
further limit the applicability of the Roche lobe model. If 
the planet has a sufficiently large magnetic field that the 
LI Lagrange point lies inside the dead zone, gas pressure 
is insufficient to open the magnetic field lines and the 
flow through the LI point is expected to be choked off. 
Also, the magnetic field may torque the gas, keeping it 
in corotation with the planet out to the Alfven radius. 
By contrast, if B 2 /8tt <C P — pv 2 at the sonic point, and 
r s,o 3* r L, magnetic stresses and tides may be ignored 
the Roche-lobe model is expected to be recovered. 

For the models of HD 209458b considered in this pa- 
per, inspection of Table Q] shows that r Sj o, Tl and the 
dead zone radius may be within factors of a few of 
each other, and the situation is more complex than the 
simplified Roche-lobe overflow model permits. 

14. SUMMARY AND DISCUSSION 

The objective of this paper was to develop a model 
for the upper atmospheres of hot Jupiters, including the 
influence of a dynamically important magnetic field. Our 
starting point (§'s 01311113 Ellwand AppendixE]) was to 
estimate field strengths for hot Jupiters, and to apply the 
theoretical model developed for MHD winds from stars 
to the case of winds escaping from the upper atmospheres 
of planets. In the process, we included strong tidal forces 
from the parent star (§'s [5j El and[7|). We computed a 
ID model of the temperature profile and ionization state 
of the atmosphere (§[8j), and constructed maps of neutral 
hydrogen column and fluid velocities to understand the 
mass loss and transmission spectra of HD 209458b (§'s|Hl 
[TOl [TTl and IT2|) . We contrast this model to the standard 
Roche- lobe overflow model (§ I13[) and verify, a posteriori, 
the validity of the MHD approximation (Appendix [B]) . 

In section [3 we discussed the application of dy- 
namo models to understand the magnetic field strength 
generated by the planet, which is currently uncon- 
strained by observation s. Using the recent results of 
iChristensen et ail ()2009[ ). which showed that the dynamo 
field increases with heat flux in the planet's core, we ar- 
gued that the large radii of hot Jupiters, and hence large 
core flux, imply that the magnetic fields of inflated hot 
Jupiters may be larger than Jupiter's field. This moti- 
vated exploring a wide range of possible magnetic field 
strengths, both smaller and larger than Jupiter's field. 

The formation of a dead zone, in which gas pressure 
is insufficient to open up magnetic field lines, was mo- 
tivated with a toy problem (§ 0) as intuition for under- 
standing the detailed structure of the hydrostatic model 
(§ [6]). The projection of the tidal force along magnetic 
field lines was used to derive the "magnetic Roche lobe 
radius" (§ [5J, outside of which gravity points outward 
along the magnetic loop. Net gravity can point out- 
ward for loops slightly larger than the distance to the 
L1-L2 Lagrange points, even in the plane perpendicular 
to the star-planet line. As a result of net outward gravity, 
the density may increase outward, as shown in Figure H 
We defined the key parameters A and e, characterizing 
the binding energy of the gas and the strength of tides, 
and their values for the observed transiting planets were 
given in Figures and [3J Many close-in planets have 
weakly bound atmospheres with A < 10, and are subject 
to strong tidal forces with e > 0.1. The magnetic field 
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strength was characterized by the plasma j3 evaluated at 
the base of the atmosphere. Solutions for the radius of 
the dead zone depend on the parameters A, /3 and e, as 
shown in Figure [5] We found that for typical param- 
eters, the dead zone extends to ~ (3 — 20)i?, implying 
that much of the volume of the magnetosphere near the 
planet is occupied by bound gas with no bulk velocity. 
Even gas outside the Roche-lobe radius can be static, if 
the dead zone is larger than the Roche-lobe radius. 

Open field lines, which are capable of supporting an 
outflow, were discussed in section [7] The momentum 
equation along field lines was used to compute the posi- 
tions of the (slow magneto)sonic points for a set of models 
using dipole geometry. Analytic solutions in the limit of 
strong and weak tides were given, which illustrated that 
inward tidal forces at the magnetic poles (for a magnetic 
dipole moment aligned with the orbital angular momen- 
tum axis) may eliminate the sonic point solutions near 
the planet. Thus, sufficiently strong tides effectively shut 
off the wind, creating a second dead zone at the poles. 
Figures [7| and [8] show solutions for sonic point radius and 
base velocity versus footpoint position. When the sonic 
point position moves far from the planet, the base veloc- 
ity becomes small, and the field lines are effectively hy- 
drostatic. Depending on e and /3, the equatorial and po- 
lar dead zones may dominate the volume near the planet. 
Lastly, we estimated the asymptotic flow speed due to 
tides in eq[40l showing that i> aS ymp "C 100 km s _1 for 
the orbital periods and stellar radii of interest. Con- 
sequently, bulk motion cannot affect the Lyman a line 
profile at Av > 100 km s _1 from line center. 

As a prelude to discussion of global models of the mag- 
netosphere, and the Lyman a transmission spectrum, 
we presented a simple spherical model of photoioniza- 
tion and thermal balance (§[8]) in order to assess the size 
of the "warm" neutral H layer. We computed the depth 
dependence of photoelectric heating in Figure |H1 showing 
that the heating drops off with pressure as a power-law, 
rather than an exponential, into the atmosphere. The 
resulting photoelectric heating, which we assumed was 
balanced by collisionally-excited Lyman a cooling, gives 
temperatures T ~ (5 — 10) x 10 3 K down to pressures 
P ~ (10 — 100) nbar. As a result, this neutral H layer 
contributes significantly to t he radius, as sh o wn in Fig- 
ure [10] As first stressed by iKoskinen et all (|2010D . the 
location of the warm H layer is key in understanding the 
large observed transit depths SF/F ~ 5 — 10%. The 
transit depth due to the layer extending upward from 
the H-H + ionization layer alone is too small to explain 
th e observations of HP 209 458b, as discussed in detail 
bv lMurrav-Clavet all (f2009h . 

Global models of the magnetosphere were constructed 
(§[9|); both to compute mass loss rates f§ [T0|) . and to con- 
struct maps of the neutral hydrogen column densities for 
a range of parameters as observed during transit (§ [TT|). 
A by-product of the warm, deep H layer is a larger mass 
loss rate than in studies w ith more shallow H layers (e.g., 
iMurrav-Clav et alJl2009h . The net mass loss rates are 
still insufficient to evaporate the planet, and are reduced 
by a factor of 3-10 due to the presence of the magnetic 
field for the parameters used. The largest columns within 
a few R of the planet occur in the dead zones, and may 
receive a contribution from H atoms outside the Roche 
lobe, but which are still bound to the planet. Hence, 



the observation of H atoms outside the Roche lobe alone 
cannot be stated as evidence for mass loss. 

The 9 global models in Table Q] were used to compute 
Lyman a transmission spectra in section 1121 We stress 
that the observational quantity most directly compara- 
ble with models of the magnetosphere is the fractional 
flux decrease between in and out of transit spectra - 
this quantity is relatively independent of ISM absorption, 
geocoronal contamination, and the background stellar 
spectrum, and is directly computable from atmosphere 
models. The compar ison between the models and data 
for HD 209458b from lBen - Jafiel (|2008f ) is shown in Fig- 
ure [TBI By variation of the base pressure of the warm H 
layer, and temperature, models can be made to bracket 
the data points, although the large error bars do not 
allow precise determination of the atmosphere's param- 
eters. Increased magnetic field is shown to increase the 
transit depth, as docs moving the planet further from the 
star. 

A comparison of the MHD wind model presented in 
this paper with the more commonly-used Roche-lobe 
overflow model was given in section 1131 It was argued 
that different regimes of accretion are possible depend- 
ing on the position of the sonic point (of an isolated 
body), the L1-L2 Lagrange points, and the size of the 
dead zone. In particular, if the LI Lagrange point is in- 
side the dead zone, gas pressure is insufficient to open up 
the magnetic field lines, and a narrow flow through LI 
is not possible. These considerations suggest that mass 
loss from hot Jupiters may be more complex than the 
simple Roche-lobe overflow model. Estimates of collision 
rates in the atmosphere (Appendix B) demonstrate the 
validity of the MHD approximation, that the e-p-H gas 
is well-coupled collisionally at the expected densities and 
temperatures in the atmosphere, and that even neutral 
H gas cannot ballistically escape the planet. 

The model presented in this paper shows that mag- 
netic fields may strongly affect theoretical estimates of 
fluid density and velocity in the upper atmosphere, and 
even the interpretation of transit depths, since neutral 
H atoms outside the Roche-lobe radius may not be es- 
caping. In future work, we hope to include additional 
physical effects, such as the interaction with the stellar 
wind, more detailed photoionization calculations includ- 
ing heavy elements, and collisional (non-MHD) effects, 
which will allow a more comprehensive physical picture 
of the upper atmospheres of hot Jupiters. 
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APPENDIX 
A. MHD WIND EQUATIONS 

In this appendix, we present the MHD equations and discuss how currents produced in the magnetosphere modify 
the field produced by the planet's core. This discussion motivates our choice of field geometry used in the global 
models. 

There is a we l l deve loped literature for axisymmctric winds from ro tating, magne tized stars. An excellent rey iew is 
given bv lSpruitl (|1996l ). Here we rely heavily on the analytic studies in iMestell (|1968f ) and iMestel fc Spruitl ()1987t ). The 
inclusion of the magnetic field can greatly affect the mass loss rate and wind speed for sufficiently fast rotation. We are 
not aware of detailed studies of wind launching from rotating magnetized bodies including the non-axisymmetric tidal 
acceleration. We postpone a numerical study of such a problem to a future investigation, here using a semi-analytic 
treatment. 

The three-dimensional MHD equations for steady isothermal flow in the frame corotating with the planet are mass 
continuity 

V-(H = 0, (Al) 

the Euler equation 

J x B 

v ■ Vv + 20 x v = -a 2 Vlnp- VC/ H , (A2) 

cp 



Ohm's law for infinite conductivity 
the induction equation 

Ampere's equation 

the isothermal equation of state 
and the no monopoles condition 



E = -vxB/c, (A3) 

V x E = —V x (v x B) = 0, (A4) 
c 

V x JB = — J, (A5) 

c 

P = pa 2 , (A6) 

V • B = 0. (A7) 

We have used constant a 2 to rewrite the pressure gradient as — Vp/p = — a 2 Vmp. The isothermal approximation is 
justified in section [5] The Coriolis and centrifugal forces appear in eq |A2l as we work in a corotating frame (see section 

To gain further insight, we rewrite eq |A2l using the vector identity v ■ Vt> = V(t> 2 /2) — v x (V x v) to obtain 



where 



Constants of the motion can be derived by dotting eq |A8l with B to eliminate the Lorentz force. We find 

B ■ VW=- (2ft + V x v) ■ (v x B) = 0, (A10) 
since the electric field vanishes in the co-rotating frame (|Spruitlll996f ). Hence W, the Bernoulli constant, is constant 



VW = v x (2fl + Vxv) + —JxB (A8) 

pc 



W=-v 2 + a 2 In p + U. (A9) 



along field lines. Another way to understand the work done on the gas is to dot eq |A8l with v 

pv ■ VW = V ■ (pvW) = --v ■ (B x J) = J E = 0. (All) 

c 

In the rotating frame, work is done on the gas by — VC/, while in the inertial frame the electromagnetic field performs 
J ■ E work on the gas (|Spruitlll996l ). 
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To understand the magnetic field structure in more detail, we take the cross product of eq |A8l with B to obtain the 
equation of trans-field force balance. Solving this equation for the component of current perpendicular to B we find 

— J± = —{J-bb-J) = \b x [VW + (2ft + V x v) x v] . (A12) 

c c v A 

Here va = B/^Airp is the Alfven speed. This equation describes the perpendicular currents that must flow in order 
to achiev e perpendicular force balance. In axisymmetry, this equation is often called the modified Grad-Shafronov 
equation (jHeinemann &z Qlberti ri978: Lo velace et al.l fl986). Perpendicular currents arise due to either vorticity in the 
flow, or variation of the Bernoulli constant from one field line to the next. In the dead zone, the fact that v = 0, and 
the further assumption that W is constant at the base, implies that J ± = in the dead zone. Parallel currents are 
determined from J ± by charge conservation, V • J = 0. 

An order of magnitude estimate for the fields SB created by volume currents Jj_, compared to the planetary 
dynamo-generated fields B p is 

B± r 4vr max(a 2 , w 2 , fire, (fir) 2 ) 

-5-~-5 — J -L 12 > ( A13 ) 



B n B c 



J A,p 



where va, p = B p j y/4np. The terms separated by commas on the right hand side of eq |A13l are estimates of the 
individual terms in eq lA12[ This estimate shows that volume currents can only significantly perturb the field out near 
the Alfven radius where v ~ fir ~ va, p - As we now discuss, at much smaller radii, of order the dead zone radius, the 
fi eld is al r eady stro ngly perturbed by curre nt sheets. 

iMestell (|1968t ) and lMestel fc Spruit! (|1987t) discussed the matching conditions between the dead and wind zones. The 
finite velocity in the wind zone acts to decrease the pressure there relative to the dead zone. Integrating the momentum 
equation across the dead zone- wind zone boundary, the total gas plus magnetic pressure must be continuous, so that 
the magnetic field strength must increase moving from the dead to the wind zone. This implies the existence of a 
current sheet separating the dead and wind zone boundaries, as shown in Figure [TJ Letting the subscripts "d" and 
"w" denote quantities just inside the dead and wind zones, respectively, total pressure continuity can be written 

p d + ir= p ™ + ir- ( A14 ) 

For identical conditions at the base, the Bernoulli equation relates the pressures as P w ~ Pd exp(— v 2 /2a 2 ). Considering 
only the dipole field from the planet, B p , and the field SB ~ 2-kK/c produced by current per u nit length, K, the fields 
in the dead and wind zones are B w = B p + SB and Bd = B p — SB. Plugging in to cq |A14l the solution for the line 
density is 

f - «-^) . (M5, 

The ratio of the field produced by the sheet current compared to that from the planet's core is then 

f^O -«-"""*') < AI6 > 

where /3d — 8irPd/B 2 is the beta for the planetary field just inside the dead zone. Inside the sonic point in the wind 
zone, v <C a and sheet currents only slightly perturb the field since SB/B ~ f3 d v 2 /8a 2 ~ (v/2va) 2 <C 1. Outside the 
sonic point, where »/a» 1, we find SB/B p ~ /3d/4, which increases outward. Hence the field configuration is expected 
to be significantly altered from the dipole outside the /3d ~ 1 point in the wind zone. As we have assumed that ^> 1 
at the sonic point, we expect the field to be altered in between the sonic and Alfven points. 

In addition to the sheet currents at the dead zone- wind zone boundary, there is a sheet current at the equator in the 
wind zone. This sheet current causes the reversal in sign of the field near the equator, approaching the split monopole 
form B cx r~ 2 e r sufficiently distant from other current sources near the planet. Since K oc B r oc 1/r 2 in the wind 
zone, and K oc (3d increases in the dead zone, we expect the maximum current to occur near the cusp in the magnetic 
field. The polar dead zone is expected to have a smooth transition from dead to wind zone, as Figure |8] shows a gradual 
transition. We expect volume currents in this transition, rather than true sheet currents. 

Based on these analytic estimates, an approximate field geometry in the wind zone is roughly dipolar inside the 
dead zone radius and roughly straight field lines outside. To go beyond this would require a detailed solution of the 
trans-field force balance for the field geometry, which is beyond the scope of the present work. 

B. MEAN FREE PATHS, ION-NEUTRAL DRIFT AND OHM'S LAW 

In this section we discuss the relative motion of the e-p-H gas as well as the magnetic field fo r the c onditions relevant 
to hot Jupiters (see Figure ITU)) . Equations and collision rates are taken from iSchunk fc Nagvl (|2004D . SN hereafter. 

To simplify the calculation, we assume all three species are isothermal with temperature T, and we work in the 
"diffusion approximation" in which inertial terms are ignored in the fluid equations of each species. Let Vj be the 
mean velocity of species j, E the electric field, and Vjk the momentum-transfer collision rate between species j and 
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TABLE 2 

Collision and gyration frequencies 



Q e = ^ = 1.8x 10 7 rad s^ 1 (B/1G) 
Qi = — = 9.6 x 10 3 rad s" 1 (B/1G) 

1 ITlpC *> ' ' 

v e n = 6.4 s" 1 (n H /10 s cm" 3 ) 
t/ Hc = 3.5 x 10- 3 s- 1 (n c /10 : 



cm 3 ) 

: 3700.0 s- 1 (n p /10 8 cm- 3 ) (10 4 K/T) 3/2 (lnA cp /10) 
: 2.0 s- 1 (n e /10 8 cm- 3 ) (10 4 K/T) 3/2 (lnA pc /10) 



i/pH = 1-2 s- 1 (n H /10 8 cm" 3 ) (T/10 4 K) 



^Hp 



1.2 s- 1 (n p /10 8 



(T/10 4 K) 



1/2 

1/2 



electron cyclotron 
proton cyclotron 
polarization 
polarization 
Coulomb 
Coulomb 
Charge exchange 
Charge exchange 



SN eq.4.88 and table 4.1 
SN eq.4.88 and table 4.1 
SN eq.4.140, 4.56 
SN eq.4.140, 4.56 
SN table 4.5 
SN table 4.5 



k. Momentum conservation implies rijrrijVjk = nkm^Vkj- We follow [Braginskii (|1965[ ) and ignore anisotropy in the 
collision frequencies, using the parallel value here for simplicity. The effective gravity be denoted g = — VJ7 and the 
pressures are Pj = njkbT. The momentum equations for e, p and H are, respectively, 

~en e + ii> e x JB^ - VP e + n e m e [g + v ep [v p - v e ) + v eH [v H - v e )\ = (Bl) 

E + ^-v p x B j - VP p + n p m p [g + u pe (v e - v p ) + v pH (v H - v p )] = (B2) 
VP H + n H m p [g + v He (v e - v H ) + v Hp (v p - v H )\ = 0. (B3) 



In order, the terms in eq IBll arc the Lorentz force, the pressure gradient, gravitational force, and collision drag force 
between e-p and e-H. We further impose charge neutrality 

n e = n p , (B4) 

and define the center of mass velocity (used throughout the paper) 

n e m e v e + n p m p v p + nHm p VH n p v p + uhVh 



m e n e + m p n p + m p nn n p + n# 



(B5) 



where the second equality is valid in the m e /m p <C 1 limit. 

The momentum transfer and cyclotron frequencies are given in Table [21 They arc shown as a function of depth in 
Figure [T71 using values of n pi uh, T and B = Bj^ eq (R/r) 3 for the hydrostatic model of the equatorial dead zone shown 
in Figure [T0l For these parameters, the gyration frequencies are larger than the e and p collision frequencies over the 
entire H + and H layers, implying motion of both e and p perpendicular to magnetic field lines is greatly restricted by 



25 



the magnetic field. Collisions with c are dominated by p well into the H layer, while H dominates collisions with p 
deeper than the H-H + transition. H atom collisions with p dominate over those from e. 

For a hydrogen atom traveling at a typical speed ch — 10 km s _1 (T/10 4 K) 1 / 2 , the mean free path against collisions 
with p is ~ ch/vhp — 10 km (10 s cm~ 3 /n p ). The proton density is sufficiently large that the mean free path is 
smaller than the scale height, ~ r 2 /XR, over the entire range shown in Figure [TU] We conclude that, due to proximity 
to the star, the high temperature and large scale height cause the density to be large enough that a fluid treatment 
is appropriate. The hot Jupiter magnetospheres discussed here are collisional, and the exobase is sufficiently distant 
from the planet to be of little practical importance. A corollary is that H atoms do not fly ballistically through the 
magnetosphere, and hence acceleration by stellar tidal gravity or radiation pressure does not cause acceleration of 
H atoms away from the planet (Lyman a radiation pr essure is only effective in a thin outer skin where Lyman a 
optical depth is less than unity dMurrav-Clav et aI1l2009D ). Rather, acceleration induces a drift velocity, which we now 
estimate. 

Ignoring the vji e term in ea lB3[ the ion-neutral drift velocity is 



vh - v f 



1 

vh p 



1 



n H m p 



(B6) 



For a simple estimate of the drift speed, we ignore the pressure gradient term, and use fiducial values g ~ 10 3 cm s~ 
and VHp ~ 1 s , giving vh — v p ~ 10 m s , and a drift time over a distance Rj of months. However, ignoring 
the pressure gradient is a poor approximation. In the H layer, hydrogen atoms provide the pressure support and 
so hydrostatic balance implies the quantity in parenthesis in cq lB6l is small. In the H + layer, in photoionization 
equilibrium, the same cancellation occurs, but for a different reason. There, both protons and electrons provide the 
pressure support, and so the proton scale height is ~ 2kbT/m p g. But in photoionization equilibrium, cq l41l implies 
riH oc n p , giving hydrogen scale height ~ kbT/m p g, so that the terms in parenthesis in eq |B6l verv nearly cancel. The 
deviations from photoionization equilibrium implies the drift velocity is proportional to a factor n p /n cq in the H + 
layer and n eq /nH in the H layer, and the drift velocity is much smaller than the naive estimate ~ g/vn p , except near 
the H-H + transition. Hence the drift time over a distance ~ R is much longer than the photoionization time of ~ hrs, 
hence photoionization equilibrium is a good approximation, as little diffusion can occur in between photoionization 
events. 

Nex t we discuss deviations from pe rfec t flux freezing. To derive Ohm's law, we follow iBraginskil (|1965| ) and solve 
eq lB3l for vh, plug the result into eq lBl| and change references frames from v e to v in the Lorentz force, with the 
result 



„ 1 „ (JxB 
E + -v x B=[ 

c \ n e ec 



P P Ph vue 
P P v h 



Ph 1 
p v H c 



— VP H xB + -+-S ff 1 + -^ 

Ph J a e \ v H 



VP e v eH V, 



vh 



Here vr = vu e 



VHp, 



{m e / n e e 2 ){v e p + v e HVHp/vH) is the conductivity, p v — m p n p , pu — m p TiH, and 



p ~ p p + Ph- The second term on the left hand side is due to induction. The terms on the right hand side are the 
Hall term, drift due to net force on the neutrals, the Ohmic term, the (small) term due to gravity on the electrons, the 
electron pressure gradient term, which gives rise to the charge separation field, and its correction due to collisions with 
neutrals. The second term on the right hand side may be put in the form of "ambipolar diffusion" , in astrophysical 
parlance, by using the total momentum equation 



1 



p H g - VP H ~ - Pp g + V(P e + P p ) J x B, 



(B8) 



yielding a term 



Bx(JxB) 

pc 2 v H 



(B9) 



on the right hand side. 

Applying cq |B7l to compute magnetic field evolution requires knowledge currents and particle densities. As argued 
in Appendix the cross-field currents are zero in the dead zone if the Bernoulli constant is uniform at the base 
of the atmosphere. While true for the simple case considered in this paper (isothermal, hydrostatic equilibrium), 
non-isothermal conditions and/or fluid motion at the base may induce perpendicular currents. 

We now discuss the relative size of terms in Ohm's law. In the H + layer, the Ohmic diffusivity is 77 = c 2 /Ana ~ 
(c 2 /47r)(m e ^ ep /n e e 2 ) ~ 10 7 cm 2 s _1 (T/10 4 K) 3 / 2 , independent of density. Assuming Ohmic decay is balanced through 
the induction term, and that the currents are of order J ~ (c/4ir)(B/r), the required (center of mass) fluid velocity 
is v ~ rj/r ~ 10 -3 cm s _1 , many orders of magnitude smaller than any characteristic velocity in the problem. 
Deep in the H layer, Figure [T7] shows that the collision rates, and hence diffusivity, may increase by an order of 
magnitude. For nonzero cross field currents, the ratio of the Hall to the Ohmic term is roughly ~ Q e /v e ~ 10 4 , where 
v e = v ep + VeH^Hp/ vh- When significant cross field currents exist, the Hall drift speed can be much larger than the 
Ohmic drift speed, but still much smaller than the gas sound speed. Lastly, if the neutral drift speed has a cross field 
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component, the second term on the right hand side may generate a fluid velocity v ~ (ph / p)(vh — v p ). This drift 
speed is much larger than the Ohmic drift speed, although it is still much smaller than the sound speed. 

We conclude that, for the ionization models discussed in this paper, the ion-neutral drift velocity and deviations 
from flux freezing in the H and H + layers are small, and single-fluid MHD is a good approximation. 
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